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FOREWORD 


This is a progress report on a research project under, 'institute for Computational and 
Applied Mechanics (ICAM)," for the period ended June 30, 1999. During this period attention 
was directed to Numerical Investigation of Hydrogen and Kerosene Combustion in Supersonic 
Air Streams’. Important results of this study are scheduled to be presented at the 9 th AIAA 
International Space Planes and Hypersonic Systems and Technologies Conference, November 1- 
4, 1999, Norfolk, VA and the 38 th AIAA Aerospace Sciences Meeting and Exhibit, Reno, NV, 
January 11-13, 2000. 

This work was supported, in part, by the Old Dominion University’s Institute for 
Scientific and Educational Technology (ISET) through NASA Langely Research Center, 
Cooperative Agreement NCC1-232. The cooperative Agreement was monitored by Dr. Samuel 
E. Massenberg, Director, Office of Education, Mail Stop 400. 


NUMERICAL INVESTIGATION OF HYDROGEN AND KEROSENE 
COMBUSTION IN SUPERSONIC AIR STREAMS 


A.A. Taha 1 , S.N. Tiwari 2 , and T.O. Mohieldin 3 

College of Engineering and Technology 
Old Dominion University, Norfolk, VA 23529 

ABSTRACT 

The effect of mixing schemes on the combusion of both gaseous hydrogen and liquid 
kerosene is investigated. Injecting pilot gaseous hydrogen parallel to the supersonic incoming air 
tends to maintain the stabilization of the main liquid kerosene, which is normally injected. Also 
the maximum kerosene equivalence ratio that can maintain stable flame can be increased by 
increasing the pilot energy level. The wedge flame holding contributes to an increased kerosene 
combustion efficiency by the generation of shock-jet interaction. 
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1. INTRODUCTION 


The scramjet engine is expected to be the most effective propulsion system, but there are 
severe conditions for ignition such as low pressure, low temperature, and high airflow speed due 
to the supersonic incoming air. Considerable fundamental research has been conducted in 
response to the increased interest in the development of scramjet propulsion systems. Many 
experimental and numerical studies of various aspects of fuel injection in the combustor are 
discussed in. 2 

The development of hypersonic airbreathing plane capable of horizontal takeoff and 
landing and acceleration to low earth orbit has attracted significant interest in the aerospace 
community for some time. It has a significant weight advantage over presently used rocket 
propulsion systems by eliminating the need to carry onboard oxygen tanks. One of the major 
challenges in developing such a plane is the design of a suitable propulsion system. Current 
airbreathing engine designs, which utilize subsonic internal flow velocities, such as turbojet or 
even RAMJET engines, would be inappropriate because of the extreme conditions seen in 
hypersonic flight. Temperature recovery from the hypersonic free stream would be in the order of 
2000 K or more, leading to engine materials difficulties and loss of usable energy due to 
dissociation of air molecules. The induced drag on the vehicle due to strong shocks at the engine 
inlets is a strong function of the flight Mach number so that it rapidly becomes impractical to use 
subsonic engines at high speeds. A supersonic combustion RAMJET or SCRAMJET has long 
been considered to be a feasible engine concept for hypersonic vehicles. 3 Figure 1.1 shows a 
schematic of a SCRAMJET engine concept where the aircraft underbody and engine designs are 
integrated together. 

♦The superscript numbers indicate references. 
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The shock waves formed by the vehicle forebody and the engine inlets serves as the 
engine compressor section while the engine exhaust and vehicle afterbody serve as the expansion 
nozzle to generate thrust from the gases burned in the combustor. The surface geometry would be 
chosen so that supersonic flow would be maintained throughout the internal passages of the 
engine avoiding large recovery of temperature and pressure from the free stream air. 4 

Figure 1 .2 shows the Airframe-Integrated supersonic combustion ramjet (scramjet) along 
with a partial cross-section through a typical modular engine. The sidewalls of the inlet continue 
the compression process, which began on the vehicle forebody. The instream struts complete the 
compression process and also provide locations for global distribution of fuel. The particular 
concept depicted here shows a combination of perpendicular and parallel fuel injection, which 
would typically take place behind a rearward-facing step and at the base of the strut, 
respectively. 5 

How these fuel-injection modes are used to control heat release is shown in Fig. 1.3. At 
lower speeds (below Mach 5, and to some extent Mach 5-7), too much heat release too early in 
the combustor will result in thermal choking and inlet unstart. Parallel fuel injection, which will 
be shown to have a slower-mixing process, is therefore used extensively to stretch out the 
combustion zone. Above Mach 7, thermal choking is much less likely to occur, and the faster- 
mixing perpendicular injection process is utilized to get faster combustion and higher 
performance. In order to quantify this mixing-controlled combustion philosophy, a considerable 
amount of research has been done on perpendicular and parallel mixing. Figures 1.4 and 1.5 
show mixing efficiency as a function of relative combustor length for each injection mode. The 
bands represent differences caused by details of the injection layout such as close or wide sparing 
of fuel injection holes, injector shape, etc. Although the curves for perpendicular and parallel 
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injection are fairly similar, they produce quite different overall results due to the fact that the 
non-dimensionalizing parameter, Xu is much larger for parallel injection. The mixing length Xi is 
also very sensitive to equivalence ratio (Fig. 1.6). This plot shows that the distance required for 
complete mixing maximizes at an equivalence ratio of unity. This is qualitatively intuitive, since 
the case where each fuel molecule must find its oxidizing molecule would clearly be more 
difficult than a case where there is an abundance of either fuel or oxidizer. The magnitudes, 
however, are very significant. Decreasing the combustor equivalence ratio to 0.8, or increasing it 
to 1.2, can cut the mixing length (and combustor length) by 30 percent or more! 

It should be noted that the data in Figs. 1.4- 1.6 are based on many measurements in real 
supersonic mixing experiments. They are, however, primarily mixing data obtained in cold-flow 
mixing situations. The impact of combustion is shown in Fig. 1.7, which shows the cold-flow 
mixing trend (cross-hatched), measured mixing efficiency (open symbols), and measured 
combustion efficiency (closed symbols) in various configurations. The abscissa is the ratio of 
combustor length to strut gap, where strut gap is the physical dimension between struts as 
depicted in Figs. 1.2 and 1.3. Note that combustion efficiency is consistently below mixing 
efficiency; this indicates chemical kinetics-related or other combustion-chemistry phenomena, 
which are potentially quite important. These issues are often related to flow-field details at the 
very local level. 


1.1 Injection Modes and Mixing 

The thermodynamic conditions at the test section entrance of the Scramjet engine are 
representative for the lower end of the hypersonic flight regime. For these vehicles it is 
anticipated that the fuel be used to cool components of the vehicle and of the engine and, thus. 
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enter the test section at high temperatures (even a supercritical conditions), resulting in fast 
evaporation and mixing upon injection. The low range of hypersonic flight regime is, therefore, 
of particular interest, since the level of heating of the fuel will be relatively low, and thus longer 
residence time will be required for the fuel to achieve complete combustion within the engine 
combustor. The mixing-combustion coupling at these conditions is of great interest. 6 

It is understood that there exists a strong interaction between mixing and combustion, 
even in the limiting cases when one or the other becomes rate controlling. Due to this coupling, 
observations from mixing in non-reacting flows are hard to extend directly to reacting flows. 
Furthermore, in many practical situations, the Damkohler number, defined as the ratio between 
the characteristic fluid residence time in the reaction zone and the chemical reaction time scale, is 
0(1). In such cases, the flowfield may become either mixing or kinetic limited in localized 
regions, however, in general, combustion takes place in thick, highly turbulent burning regions. 
The mean and time varying parameters responsible for mixing are affected by heat release, 
which, in tum, affect the chemical kinetics; thus, mixing and combustion become closely 
coupled. 7 

Different methods of liquid fuel injection have been investigated. Parallel flow mix 
through the growth of the shear layer generated at the interface of the different components. It is 
observed that compressibility plays a major rule in reducing the growth of the shear layer, thus 
requiring longer fluid residence time to insure good mixing in practical devices. As the 
compressibility increases, the spreading rate of the shear layer was observed to drop to about 14 
of the compressible shear layer growth for the same velocity and density ratios of the two mixing 
constituents. The convective Mach number (defined as a measure of the speed of the turbulent 
structures in the shear layer relative to the free stream) increases the stability of the large-scale 


4 


turbulent structures which are responsible for the turbulent shear layer growth. This is of 
particular significance in high-speed flows, where parallel injection is the preferred mode due to 
the lower total pressure relative to the other fuel injection and mixing options. An increase in the 
stability of the large vortical structures in the shear layers due to compressibility is observed even 
in situation where the turbulence is intentionally intensified. Growth of the shear layer was 
forced by a pulsating shock impinging on the shear layer in parallel compressible flows in high 
speed (i.e., Mach 3 and 5, respectively) as a destabilizing mechanism leading to the growth of the 
shear layer and accelerated mixing. Only a negligible effect was obtained both in the near and far 
field. However, when a shock of the same strength impinged on the boundary layer of the two 
flows prior to their coalescence, an increase in the shear layer growth (14 to 26 percent) was 
experienced. Thus an increase in the turbulence in the boundary layer upstream of the mixing 
station increased, causing an increase of as much as 26 percent in the shear layer growth. This 
indicated the significant role played by the small-scale turbulence in the mixing process. 
Introduction of axial vorticity was shown to improve mixing, which was also attributed to the 
small-scale turbulence effect on the growth of the shear layer and on the enhancement of 
turbulence in the shear layer. Mixing at a molecular level is, in particular, important in 
combustion applications when stoichiometric mixing is a prerequisite to initiate chemical 
reactions. 6 

The effect of small-scale turbulence on the mixing process in high-speed flows requires 
more detailed investigation. In particular, mixing at a molecular level, which indicates the ability 
to initiate chemical reactions, is of a great interest. Concentration distribution is not sufficient to 
provide this information. Although it is recognized that a reacting flowfield is substantially 
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modified once exothermic reactions take place, mixing of non-reacting flows needs further 
investigation to provide further insight of the characteristics of the micro-mixing processes. 6 

Studying the flow field behavior of the transverse injection of the fuel is one of the recent 
demanding research trends in the combustion field of the scramjet engine. This technique is used 
to increase the fuel-air mixing in order to achieve the required heat release pattern with the short 
combustion residence time associated with the high Mach number condition. 2 

Transverse or inclined injection introduces different mixing mechanism. At a macro- 
scale, axial vorticity is produced at the edges and above the jet as the main flow surrounds the jet 
in an axial and transverse direction. These vortices break the plume and entrain the injectant into 
the core flow. The level of mixing obtained by this strong momentum exchange is traditionally 
quantified by the degree of penetration and concentration decay of the jet injected into the main 
flow, usually air. 6 Studies performed by Mays et al. 8 and Wood et al. 9 addressed the effect of 
small-scale turbulence on transverse injection and mixing. It was suggested that the onset and 
divergence of an instability in the jet is ultimately responsible for the plume fracture. This 
instability propagates upstream and produces additional turbulence in the core flow. Even in 
supersonic flows, an upstream interaction exists, as the instability at the jet boundary induces 
oscillations of the bow shock generated in front of the jet via the separated region in front of the 
jet. In tum, these shock oscillations contribute to increase the vorticity in the flow. 6 

When heat release effects are included, significant modifications of the turbulent 
structures occur. It has been indicated that the heat released from chemical reactions both 
generate and suppress turbulence. It is shown that the turbulent dilation and viscous dissipation 
in flames reduce small-scale turbulence, while the shear gradient effects, mainly of pressure and 
density generate new large-scale structures. 6 
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Figure 1.8 shows a schematic of the region surrounding a single perpendicular injector. 
Because the mainstream flow is supersonic, there is a bow shock off of the underexpanded fuel 
jet and internal wave structure associated with the jet expansion. The boundary layer ahead of the 
jet separates, and there is a recirculation zone downstream of the injector as well. The details of 
these local zones are quite important to the ignition, flameholding, and combustion processes and 
therefore important to real combustor design. 5 

Mohieldin and Tiwari 10 studied numerically the advantage of using the tandem injection 
on the regular single step tangential injection regarding the mixing aspect in the flow field. 
Matsuo and Mizomoto 1 studied the flow structure of supersonic flow of Mach 2.5 past a 
backward- facing step with perpendicular sonic injection from the bottom wall. They examined 
the effect of both the step height and dynamic pressure ratio of the jet flow to the main flow 
using the perfect gas condition in two-dimensional space. The effect of H 2 -transverse injection 
on the characteristics of the flow field was investigated numerically by Berman et al. 11 using 
both adiabatic and constant temperature wall boundary conditions. 

Schetz et al. 12 presented a comprehensive review of the mixing of transverse jets and 
wall jets in supersonic flow. While the stream wise injection has the advantage of adding to the 
thrust component of the engine, many of the approaches utilized to improve the wall injection 
have been used successfully to enhance the fuel-air mixing. 2 

Using the rearward-facing step in the supersonic combustors and injecting the fuel 
perpendicular or at least inclined on the main air flow have the advantage of holding the flame 
because of the existence of the separated and recirculating regions in the flow field as depicted in 
Fig. 1.8. This can help solving, in part, the complication of holding the flame in supersonic 
engines due to the fact that the flow in such engines is not premixed which means that if bluff 
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bodies are to be used as flame holders they have to be large. This is highly undesirable in 
supersonic combustors due to very high drag. Therefore, this problem of flame holding could be 
solved using the rearward-facing steps in the combustors and/or the perpendicular fuel injection. 
This issue emphasizes the importance of flow field details around the fuel injectors. * 

An important feature of the normal injection flow field is the gross penetration of the jet 
into the flow. Baranovsky and Schetz 13 introduced the injection angle as a separate parameter in 
the correlation, which was set to evaluate the penetration of liquid jets. It was found that the 
upstream injection (injection angle is greater than 90° with respect to the incoming air flow 
direction) can be used as a method to increase penetration and residence time of a liquid fuel into 
a scramjet combustion chamber. It was also found that the maximum penetration takes place at 
injection angle (0) =135 degrees. 


1.2 Ramp Injection and Mixing 

Enhanced fuel/air mixing schemes are required in SCRAMJET engine designs because of 
the short residence time (on the order of ms) of the gases in the supersonic flowfield of the 
combustion chamber. Uniform and rapid mixing allows for shorter combustor lengths, higher 
combustion efficiencies and better engine stability. Various mixing enhancement schemes have 
been proposed including contoured injection orifices, flowfield and shock oscillations, baroclinic 
vortex generation, and wall-mounted vortex generators. 4 

Figure 1.9 shows schematic representations of the swept and unswept ramp fuel injector 
configurations investigated in the study of Donohue. 4 The injector is located on one wall at the 
inlet of the combustor test section. Fuel is injected from the base of the ramp in the downstream 
direction, normally with some component upward to increase the penetration of the fuel into the 
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freestream air. The ramp sidewalls can be swept back to improve the mixing performance of the 
ramp. A single injector is used for the studies performed by Donohue 4 although normally an 
array of injectors would be used on several walls and possibly on struts crossing the flowfield in 
order to fuel the entire cross section of incoming air. The single ramp geometry used can be 
thought of as one element of an array with symmetry conditions at the top and side walls. More 
complex injection configurations incorporating transverse injection for better penetration, steps 
or blunt bodies for flame holding, and wall expansion or compressions would probably be used 
in conjunction with the ramp injectors in a final design to allow for design optimization at a 
variety of operating conditions. 

A diagram of the important flowfield features in the ramp injector flowfields is shown in 
Fig. 1.10. The ramp deflects the supersonic free stream, compressing the air along its top surface 
and generating a shock wave, which reflects off the opposite wall of the wind tunnel. The higher 
pressure air on top of the ramp spills over into the relatively lower pressure region to the sides of 
the ramp generating a pair of counter-rotating streamwise vortices. Fuel injected from the base 
region of the ramp is mixed into the freestream air by the ramp-generated vortices, which are 
convected downstream along the sides of the fuel jet. The side walls of the swept ramp are 
angled back to expand the flow at the sides of the ramp, decreasing the pressure in this region 
and increasing the strength of the vortices formed. The ramp base area provides a subsonic, high 
temperature region in the flowfield important for flame holding in the reacting case. 

There are a number of studies specifically related to the 3-D mixing flowfield of ramp 
injectors. Hartfield et al. 14 presented measurements of injectant mole fraction in the flowfield of 
a swept ramp with air injection using the same technique used by Donohue. 4 The results 
illustrated the strong effect that ramp generated vortices have on the mixing process. Waitz et 
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al . 15 presented both numerical simulations and experimental measurements of total pressure and 
fuel concentration for helium injection behind an unswept ramp which added the effect of 
baroclinic vortex generation when a low density fuel is used. Riggins et al. 16 presented 
numerical simulations of the mixing flowfields of both swept and unswept ramps with hydrogen 
injection into high temperature air and went on to calculate the reacting flowfield as well using 
the SPARK code. These results showed increased penetration and spreading of the jet plume in 
the reaction case. Daso and Gross 17 computed the flowfield behind swept and unswept ramps 
with air injection using the USA code and illustrated how side sweep of the ramp significantly 
increased the mixing downstream. A numerical study designed to investigate the relative 
importance of ramp generated vortices vs. baroclinic vorticity on supersonic mixing behind 
swept ramps was performed by Donohue et al . 18 

1.3 Autoignition 

McClinton 19 showed that combination of the step-base and the perpendicular injection 
downstream of the step was effective in enhancing the ignition ability. Huber et al. 20 suggested 
possible autoignition source in the supersonic combustors and proposed a simple model to 
predict ignition limits in typical combustors. These candidates are, the base region behind steps 
and struts, recirculation regions upstream of the fuel jets injected perpendicular to the airflow, 
and the stagnant region behind the bow shocks upstream of the fuel jets. 

Tomioka et al. investigated the flow field at pre-ignition phase in a supersonic 
combustor with perpendicular injections behind a backward-facing step to understand the 
mechanism of the ignition enhancement observed in the ignition tests with the same combustor. 
It was concluded that the interaction of the separation region upstream of the fuel jet and the two- 
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dimensional step base-recirculating region, namely the merging of these two regions, caused 
enlargement of the ignition region and the enhancement of the ignition ability. Ignition 
parameters were compared for the cases with and without the step. As results, the interacted 
region in case with the step was found to be more preferable ignition source due to its enlarged 
size. 


1.4 Liquid Hydrocarbon Fuels 

1.4.1 The Problem of Cooling in High-Speed Vehicles and The Role of Fuel 

As air breathing-engine-propelled vehicle speeds increase, thermal problems multiply 
because of the effect of stagnation temperature. While total cooling needs increase, the most 
critical regions are the leading edges and the engines. Although thermal effects can be somewhat 
accommodated by improved materials and passive cooling, sustained hypersonic flight in the 
atmosphere requires a substantial heat sink. Compared to a mechanical refrigeration system or a 
non-combustible coolant, the fuel is the best source of cooling. 22 

In the high Mach number region of the flight envelope (M>4), aerodynamic heating is too 
great for conventional structural materials to survive without active cooling. The fuel not only 
must have a good heat-of-combustion but also must provide the necessary heat sink for the 
cooling system. 23 

Higher aircraft speeds also have a direct impact on the operating environment a jet fuel 
will encounter. The higher speeds mean higher air stagnation temperatures, which increases the 
aircraft cooling requirements and prevents the use of air as a coolant. Thus, increasing engine 
thrust-to-weight ratio and aircraft speed result in large heat loads that must be managed with the 
main coolant available on the aircraft-the fuel. 24 


11 


The fuel in modem military aircraft is the primary coolant for on-board heat sources. It is 
used to cool aircraft components such as the engine lubricant oil, hydraulic fluid, environmental 
control system (ECS), avionics and electrical systems, and at high Mach numbers, the air frame. 
Figure 1.11 illustrates the heat loads from the different aircraft components. As aircraft and 
engine technology have advanced from the F-4 to tomorrow’s advanced fighter, the heat sink 
requirements of the fuel have increased significantly. The engine is the main heat source for the 
aircraft at Mach 3 and below, but the ECS is becoming an increasingly important part of the heat 
management problem. Fuels used in high-speed aircraft will have to absorb large amounts of 
excess heat for aircraft thermal management purposes. This will result in the fuels being heated 
to supercritical conditions. 24 

There are some key technical problems associated with flight in the intermediate Mach 3 
to 7 range that still need to be addressed. Nearly all of those technical challenges have to do with 
the handling of the excess heat that is produced either by the propulsion system itself, or by the 
aerodynamic heating of the airplane surface. Mechanical systems, electric systems, structural 
elements, passengers, and crew all require cooling. Preliminary engineering estimates suggest 
that these cooling tasks cannot be performed using conventional aircraft design approaches. 
Traditionally, this "thermal management" challenges has been met by using the onboard fuel to 
absorb excess heat. Fuel is circulated through hot sections of the aircraft, usually back into the 
tanks, and finally to the engine where it is burned. This "recuperative" approach has the benefit 
that thermal energy dumped into the fuel is then eventually recovered as additional thrust. 25 
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1 . 4.2 Supercritical Conditions and Their Effects 

Definition of supercriti cal conditions; The critical point of a fluid marks the terminus 
of the vapor liquid coexistence curve. A fluid is said to be "Supercritical" when its temperature 
and pressure exceed the temperature and pressure at the critical point. 26 

Effect of supercriti cal conditions on fuel properties: If the fuel temperature becomes 
too high, the fuel will chemically decompose to form gums and solids that cause sticking of the 
fuel control valves and fouling of the fuel nozzles and heat exchanger. 

It is clear that the future fuels must operate effectively at higher temperatures otherwise 
there will be serious limitations to aircraft performance. It is estimated that 900 F is the highest 
operating temperature that can be expected for an advanced hydrocarbon fuel without thermal 
decomposition (JP-900). At higher Mach numbers (>5), more cooling is required than can be 
obtained with JP-900 fuel. Such aircrafts will require an "endothermic" fuel, where additional 
heat is absorbed by an endothermic reaction of the fuel. 

For speeds above Mach 8, the heat loads are so large that "cryogenic" fuels will be 
required. The physical and chemical properties of supercritical fluids (SCF) are expected to be 
dramatically different from gases and liquids. It has been found that reaction rates and 
mechanisms in supercritical fluids can be quite different from those in subcritical fluids. 
Physical data for larger hydrocarbons (that might be expected to be part of an advanced fuel) is 
sparse. 24 

For very high temperature test where cracking -decomposition of the fuel into lighter 
components, which will be discussed later- has occurred, the formation of light gases may 
significantly modify the fuel properties. It appears that n-dodecane (C 12 H 26 ) properties are similar 
to those of JP-5 (Cn.9H22.2X Jet A (C116H22), and JP-7 (C12.1H24.4) and could be used in heat 
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transfer calculations. The presence of cracking products can significantly modify the fuel’s 
physical properties. At pressures near (but above) the critical pressure the physical properties (C p 
especially) are very sensitive to temperature. 24 

Thermal stability o f the fuel: The understanding of the thermal stability of jet fuels at 
supercritical conditions is almost nonexistent. There is some deposition data for fuels that 
includes supercritical conditions. As illustrated in Fig. 1.12, typically the deposition in "heat- 
tube" tests is seen to peak at about 700 °F (a range which includes the critical temperature for 
typical kerosene fuels), drop off as temperature is increased further, then begins to increase again. 
The deposition peak is attributed either to a change in mechanism or to changes in fuel properties 
as the fuel is heated to supercritical temperatures. Most hydrocarbons tested had similar critical 
temperatures. Recent results at Wright Laboratory indicate that the behavior in Fig. 1.12 is due 
more to the temperature dependence of chemical reactions than to the critical properties of the 
fuel. The fact that the supercritical deposition mechanism is unclear is not surprising since the 
mechanism controlling solids deposition in subcritical fuels is not totally known. 24 

Supercritical fuel injection: It might be expected that the injection and mixing of a 
supercritical fuel might resemble gas injection more closely than liquid injection. It was reported 
that compressibility can affect atomization in an analytical model. 27 

An additional complication is that the heat added to the fuel (and thus the fuel’s 
temperature) varies throughout an aircraft’s mission. Thus the fuel might be injected as a liquid 
during the acceleration portion of a flight and as a supercritical fluid during cruise and/or decent. 
It has been found that two-phase fuel injection can dramatically affect spray behavior. 28 
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1.4.3 Endothermic Fuels, The Promising Fuel for Hypersonic Combustion Systems 


Definitions and differences between cryogenic and endothermic systems: The 

"Cryogenic" fuels are the fuels that have a negative, less than zero, "Boiling Point". Therefore, in 
order to use these fuels as liquids, they must be cooled beyond the boiling point and storage 
environments must satisfy this condition. Examples for the cryogenic fuels and their 
corresponding boiling point temperature are the liquid hydrogen H 2 (-423 F) and liquid methane 
CH4 (-259 F). 

The "endothermic" fuels are the fuels that extract heat while forming from their initial 
components. Therefore, the decomposition of these fuels into the initial components requires 
energy for its accomplishment and by that a quantity of heat can be absorbed from the 
surrounding environment causing the required thermal loading management of the vehicle. An 
example for the endothermic fuels is the methylcyclohexane MCH (C 6 H 1 CH 3 ). The boiling point 
of MCH is 213 °F which means that it is a liquid in the normal ambient conditions. 

Cryogenic fuels: advantages and disadvantages: The cryogenic fuels have high energy 
contents and thermal tolerance. Figure 1.13 shows that hydrogen has almost three times higher 
energy per unit mass than other fuels. 29 

The temperature limit of a fuel is a key consideration for its selection as a high speed 
transport fuel. Figure 1.14 shows that hydrogen has the highest temperature limit over other 
fuels. On the other hand. Fig. 1.15 indicates that cryogenic fuels have relatively poor energy 
contents per unit volume which is an undesirable property that is magnified by their storage 
requirements . 29 

Cryogenic fuels must be viewed with extreme caution because of the severe constraints 
that they place on the design of the aircraft and because of the high cost of transporting, storing, 
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and delivering them. Hydrogen or methane handling facilities at just a few U.S. airports would 
be very costly to install and maintain, and of course these fuels raise serious safety questions. 
From the military standpoint, cryogenic fuels appear impractical because, again, military 
operations must not be tendered by a need for costly and exotic infrastructures, and by all means 
not by an uncertain and vulnerable supply of mobility fuels. Nevertheless, serious consideration 
is being given to such option and it is possible that successful large-scale applications using 
hydrogen or methane might eventually be developed. 8 Figure 1.16 shows that cryogenic fuels 
will require totally new fuel delivery, distribution and storage facilities, as well as extensive gate 
and terminal modifications. 29 Figure 1.17 shows that there will always be losses for cryogenic 
liquid during refill and tank to tank transfer processes. 29 

Endothermic fuels; the phenomenon and advantages: An innovative approach based 
on more conventional types of fuels takes advantage of the fact that some hydrocarbons, when 
they degrade at high temperatures, do so cleanly, i.e., without forming carbonaceous surface 
deposits. Instead, they break down to smaller hydrocarbons that are then circulated to the 
combustor. The chemical breakdown itself absorbs heat, thus providing another useful thermal 
management process. A number of these "endothermic" fuels are known, and some are under 
serious study. All of the known endothermic fuel reactions require a catalyst to make them 
proceed and absorb extra heat at moderate temperatures. Hence, variants such as the catalytic 
heat exchanger are being considered in future aircraft concepts. While the endothermic fuel 
technique introduces a new element of complexity into aircraft design (and particularly new 
maintenance requirements and reliability concerns), it has obvious advantages over the use of 
cryogenic fuels. 25 There are indications that, free of certain environmental and chemical 
influences, most hydrocarbons tend to have good high-temperature tolerance. 25 
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Cryogenic fuels can contribute only sensible and latent heat, whereas certain hydrocarbon 
fuels can in addition provide cooling through endothennic reactions. Hydrocarbons can undergo 
both thermal and catalytic reactions. Theoretically, the total heat-sink capacities for hydrocarbon 
fuels range from about 50% to 112% of the cooling capacity of hydrogen (on a "%" heat of 
combustion basis) with laboratory proven capabilities up to about 85 %. 22 

The endothermic fuels have the following technical, commercial, and safety advantages; 

• Continued availability 

• Relatively inexpensive 

• Condensed phase under standard temperature and pressure (STP) 

• Higher latent heat of vaporization 

• Safe handling and storage 

• High energy content per unit volume 

• No ignition limitations under subsonic combustion conditions 

• Non-cryogenic; ** insulation not required ** vapor recovery not required 

• Conventional Army/Commercial (A/C) fuel system 

• A/C turn-around time not governed by fuel 

• Conventional fuel handling/logistics; ** off-site production 

• Launch side hardening easier 

• Available enabling technologies 

• Potentially smaller vehicles 

The Cooling capacity of normal and endothermic hydrocarbon fuels: Most 
hydrocarbons can provide a maximum of only 700-800 Btu/lb enthalpy on being heated from 
ambient to about 1000 °F. Figure 1.18 shows the heat sink in terms of Btu/lb fuel burned as a 
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function of Mach number. The indication being that about 2000 Btu/lb will be required at a speed 
of Mach 10. It is obvious that flight speeds greater than about Mach 5 will require additional 
cooling capability. This region generally known as the hypersonic region, is where endothermic 
reactions will have to be utilized in order to employ the higher hydrocarbon fuels successfully. 22 

Types of he at-sink reactions (endothermic reactions); 

• Thermal reactions (Non-catalytic reactions) 

• Catalytic reactions: 

a) Dehydrogeneration reactions 

b) Dehydrocyclization reactions 

c) Depolymerization reactions 22 

Thermal reactions: The decomposition of hydrocarbons requires energy for its 
accomplishment, i.e., it is a possible heat-sink reaction. The effect of product distribution for the 
thermal decomposition of some hydrocarbons is shown in Table 1.1. At higher conversions, the 

actual reaction is essentially thermally neutral due to increased formation of methane (which is 
exothermic). 22 

Catalytic reactions: A catalyst may serve either or both of two purposes; to increase the 
rate of a particular chemical reaction or to cause one of several possible reactions to occur 
selectively. It is very useful to increase the selectivity of the endothermic reaction. For example, a 
catalyst that would selectively enhance the rate of certain endothermic cracking or 
dehydrogeneration reactions of propane to eliminate the formation of methane could increase the 
heat sink by as much as 300% over the thermal reaction. 22 
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A Dehydrogeneration reactions 

The most practical way to do this is to have the reaction go selectively to a 
product having a resonance-stabilized structure. This not only will produce a more 
reliable heat sink but also will influence the temperature regime of 
thermodynamic equilibrium in favor of the desired products. 

The ultimate fuel may be a blend of hydrocarbons in order to achieve 
optimization in such properties of freezing point, volatility, lubricity, heat sink, 
thermal stability, etc. For example, addition of methyl groups to he basic ring 
structures can alter these properties, but the more methyl groups that are added the 
lower is the available heat sink. 22 
A Dehydrocyclization reactions 

In principle, any hydrocarbon structure with sufficient carbon atoms in it could 
be cyclized to from C 5 or C 6 rings and then dehydrogenated to give a high heat- 
sink. For example, normal octane can be cyclized to ethylcyclohexane and then 
dehydrogenated to styrene with a heat sink of about 400 Btu/lb. However, the 
cyclization step on any catalyst investigated so far is relatively slow and not too 
selective, which means that the theoretical heat sink cannot be realized. It appears 
that the most promising utilization of this type of reaction will come in connection 
with already-formed cyclic structure. 22 
A Depolymerization reactions 

Since the chemical reaction of a polymerization is exothermic, it follows that 
the reverse reaction of depolymerization is endothermic. This reaction provides a 
limited heat sink . 22 
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Comparison of thermal and catalytic reactions; In general, the catalytic reactions are 


more favored rather than thermal reactions. 

# Concerning the degree of conversion : 

• The heat sink provided by the catalytic reaction is a direct function of the degree of 
conversion. The greater the degree of conversion, the greater the heat sink. 

• The thermal reaction has a maximum heat sink of about 300 Btu’s, which occurs at about 
60% conversion, and thereafter the heat sink actually decreases because this reaction is 
not selective and there are competing endothermic and exothermic reactions. 

# Concerning the rate of reaction (K): 

• The catalytic reaction has a high intrinsic rate and a low energy of activation and, hence, a 
high rate of reaction over a broad range of temperature. 

• The thermal reaction on the other hand, proceeds very slowly at low temperatures. 
Although, because of the high energy of activation it responds strongly to the effect of 
temperature it does not reach comparable rates to the catalytic reaction until about 1200 
°F. This means that the catalytic reaction has a much broader range of temperature over 
which it can be used successfully for cooling. 

# Since the rate of reaction is primarily a function of the properties of catalyst, it allows for the 

possibility of increasing the rate by devising better catalyst. This is not possible of course with 
thermal reactions. 22 

How can endothermic fuels accommodate the varyine heat load? For endothermic 
fuels, this varying heat load may also be translated to variations in fuel phase. For example, the 
tow types of endothermic fuel reactions create mixtures of small and large hydrocarbons, these 
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reaction types are: dehydrogenation such as methylcyclohexane => toluene + H 2 , and cracking, 
such as C„H 2 n +2 => light gases such as C 2 H 4 . 

These fuel reactions are used to absorb aircraft heat by passing some portion of the fuel 
through a reaction. Thus, during early portions of a mission, the fuel would probably by-pass the 
reactor and being injected as a liquid. For high heat load parts of the mission, the fuel might pass 
entirely through the reactor, emerging as a supercritical fluid. For intermediate heat load parts of 
a flight, part of the fuel might by-pass the reactor. The fuel to be burned could then consist of 
part unreacted fuel (liquid) and part hot, reacted fuel (SCF). 

The fuel could be injected into the combustor through two separate sets of nozzles 
(one for liquid and one for SCF), or the fuel could be re-mixed and injected. 24 

Cooling application possibilities: It is possible that turbine engines would be used up to 
a speed of about Mach 4, a subsonic combustion ramjet at speed up to Mach 7, with a supersonic 
combustion ramjet being used at speeds above that. 22 

Three different ways of transferring the heat from the engine to the reactor exchanger 
have been considered. These three ways are illustrated in Fig. 1.19, which represents ramjet 
engines traveling from right to left; the procedure is described below. 

• Direct regeneration cooling between the combustion space and the fuel tubes 
(containing catalyst) forming the engine. 

• Secondary loop design in which a circulating fluid is employed to transfer the heat 
from the engine to a reactor-exchanger containing the catalyst. 

• A system in which some of the inlet air (say 5-10%) is bled through a reactor- 
exchange system and then employed to cool the walls of the engine by transpiration or 
film cooling. 22 
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A Proposed cooling system ; This cooling system design was driven by Mach 5 cruise 
conditions. Operating conditions for Mach 5 are shown in Fig. 1.20. The maximum pressure is 
424 PSIA (2.93 MPA) at the inlet to the cooling panels, and the maximum temperature is 838 
°F (448 ° C ) at the outlet of the cooling panels. The cooling panels are arranged in senes with the 
coolant flowing forward through the airframe nozzle, then downward to cool the turbofan inlet 
close-off door seals and back through the ramjet. The 77 Ibm/s (35 kg/s) flow splits after the 
last cooling panel outlet with 40 lbm/s (18 kg/s) flowing through the heat exchanger/reactor and 
the rest flowing to the pre-heater. The flow is re-mixed at the inlet to the expansion tank. The 
pump requires 143 HP at 242 PSIA (1.67 MPA) and 80% efficiency. A separate power unit 
using LOX and MCH drive these pumps. The use of cryogenic liquids in large quantities is not 
allowed on aircraft carriers; however, small amounts of LOX are used for pilot breathing at 
high altitudes. It is estimated that the quantity of LOX required for the APU would be 
comparable to that required for the pilot. 23 

Syltherm 800 has a high coefficient of thermal expansion, making an expansion tank 
necessary. Between room temperature and 838 °F, Syltherm expands 91% in volume. According 
to the manufacturer, this tank should be designed to be 1/4 full of coolant when cold and 3/4 full 
when hot. 5 The remaining volume is filled with a nitrogen gas blanket. The tank is located at the 
lowest pressure point in the system, and the constant flow of coolant liquid through the tank 
ensures that vapors will be isolated in the tank. The physical arrangement of the components is 
shown in the scale inboard profile in Fig. 1.21. 23 

A unique feature of this cooling system is that it is self-compensating, i.e., an increase in 
heat load will lead to an increase in the temperature of the coolant entering the heat rejection 
devices. The increase in temperature will lead to an increase in the efficiency of the heat 
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exchanger/reactor, increasing the temperature drop in the coolant. Similarly the preheater will 
also transfer more heat to the fuel. The effect of the increase in heat load on the temperature of 
the coolant will have been diminished and a new equilibrium will be established. In this off- 
design case, the coolant and catalyst may be degraded, but the system will operate until the fuel 
temperature reaches 925 °F (496 °C); this is the limiting condition for both the catalyst and 
coolant. Syltherm can operate for five days at a temperature of 932 °F (500 °C), but the pressure 
would have to be increased to avoid boiling. Note that the heat capacity of the coolant is so great 
that a 15% increase in heat load will only cause an increase of 6 °F (3.3 °C) in the coolant. 23 

1.5 Piloted Energy for Liquid Hydrocarbon Combustion 

Liquid hydrocarbons are attractive candidates for the low range of the hypersonic fight 
regime due to their high volumetric energy content and the relative simplicity of the operational 
logistics. Although the energy/mass density of liquid hydrocarbon is lower than that of hydrogen, 
some of the mass increase wuold be recovered by a small and ligher structure of the vehicle. 
Liquid hydrocarbon fuels require substantial residence time to achieve vaporization and complete 
exothermic reactions, and such time is unlikly to be available in reasonable-sized supersonic 
combustors. The chemical kinetics of hydrocarbons are slow in comparison with gaseous fuels, 
such as hydrogen. For many realistic operating conditions, exothermic reactions cannot be 
achieved within the available residence time. Use of a pilot flame with dast kinetics (i.e., gaseous 
hydrogen) can provide locally, at the liquid injection location, the conditions necessary to 
accelerate the hydrocarbon reactions rates and reduce the Damkohler number, i.e., high 
temperature and low fluid velocities. For such a system it is necessary to verify: (a) the level of 
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piloted energy required to ingite and maiantain stable comustion of the liquid fuel, and (b) the 
interaction between the pilot flame and the heat sink represented by the injected hydrocarbon. 6 

Bonghi et al. injected liquid toluene at 5H down stream of a step of height H in a Mach 
2.5 flow at 300 K and 1000 K using a parallel hydrogen-pilot flame. Toluene (C 6 H 5 CH 3 ) is on of 
the components of catalytically cracker methylcyclohexane (MCH), a candidate fuel for 
supersonic combustion applications. Segal and Young 6 using the same previous configuration 
injected the toluene from a reservoir at 300 K, and thus, absorbed energy from the surroundings 
to vaporize and further heat to the local temperature. The large recirculation region formed by 
hydrogen-pilot combustion extends beyond the liquid injection station, 5H, and thus a substantial 
amount of the injected liquid reaches the region of stabilization of the pilot flame. When the 
amount of liquid increases above a certain quantity, it quenches the pilot flame. The amount of 
liquid fluid that can be injected without inducing quenching was found to be dependent strongly 
on the equivalence ratio of the pilot flame and, in smaller proportion, on the air stagnation 
temperature. It was found that as the pilot energy level increases the amount of toluene that can 
be injected (while maintaining a stable flame) can also be increased. Beyond a certain value, 
however, the ratio of energy levels decreases again, as more liquid is injected and the amount that 
reaches the region of stabilization of the pilot flame increases, inducing quenching. This effect is 
enhanced by the effect of heat release on the structure of the flow field. A larger deposition of 
heat in the combustion region caused a pressure increase followed by an enlargement of the 
recirculation region and a reduction in the local velocity, which favors upstream interaction, 
further facilitating the arrival of liquid in the recirculation region and quenching of the pilot 
flame. The increase in the stable flame margin with the air stagnation temperature, T 0 , was 
evident. 
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1.6 Importance of Modeling and Simulation 

The supersonic combustion ramjet (scramjet) is expected to be the most effective propulsion 
system for the Single Stage To Orbit (SSTO) transportation vehicles and hypersonic 
transportation vehicles of the next Generation. Many studies on scramjet components such as 
inlet, combustor and nozzle have been carried out to obtain a better understanding of the 
performances of the individual components. 31 ' 34 However, it is expected that intensive 
interactions among these components will occur in real engines. For example, a rise in pressure 
due to combustion causes separation which propagates upstream into the inlet and changes the 
inlet back pressure, which may result in unstart condition in the inlet. In the case of interactions, 
the total performance of the engine cannot be evaluated based on the linear combination of the 
obtained performance of the individual components. Thus, tests of the models of the whole 
engine are necessary to elucidate the interactions among these components and the overall 
performance of the whole engine. 35 Testing of the models of whole engine requires rather big 
and expensive wind-tunnel facilities. Thus, only a limited number of results on whole engine 
tests have been reported. 36,37 

Multiphase combustion is such a complex phenomenon that it warrants focused research 
using modeling and numerical simulations. Though physical experiments are the ultimate test to 
study the performance of combustion/propulsion systems, they are often extremely expensive and 
complicated, and at times are not even possible. This is due to the hostile environments, 
complicated interactions, and couplings. On the one hand, these couplings make it impossible to 
study the effect of one parameter at a time in a physical experiment of this nature. On the other 
hand, numerical experiments open an avenue for isolating and understanding the influence of 
different parameters on complex combustion situations. Furthermore, these can serve as means of 
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evaluating scaling laws, which would be extremely difficult and expensive to determine from 
physical experiments alone. Modeling also provides a potential design tool for future combustion 
and propulsion systems. However, this is not a simple or straightforward extension of 
computational fluid dynamics (CFD) by including a few more terms in the governing equations. 38 

Numerical simulation is time consuming even with ultrafast modem computers. The fuels 
of choice, e.g., hydrocarbons or metal slurries, will have a large number of kinetic steps and 
species involved in their combustion. The computational time increases as the square of the 
number of species, and calls for reduced chemistry based on sound logic. The presence of solid 
particles, such as soot or unbumed fuel droplets, and their transport in the turbulent flow field 
complicates the situation even further, and requires extremely small grid sizes. Ultimately, the 
results of these complex numerical simulations should be made available in the form of user- 
friendly models for the designer. 

Strong interaction of the computational combustion community with experimentalists is 
needed in order to assure meaningful model validations, and to generate more reliable and 
advanced models to aid in the design of complex combustion/propulsion systems of the future. 

Computational fluid dynamics (CFD) simulations have become invaluable in the design 
process due to the high costs and long turnaround times encountered in experimental programs 
and due to the inability to reproduce all desired flowfield conditions in present laboratory 
facilities. Fundamental numerical techniques for the calculation of compressible flowfields have 
been available for some time, but the computational resources required for making accurate 3-D 
simulations of complex flowfields have become more generally available for design purposes 
only relatively recently. This is, for the most part, due to the many dramatic improvements that 
have occurred in the performance of computer hardware, but is also due to progress in the 
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development of efficient and accurate numerical techniques and algorithms. Advances in areas 
such as grid generation, development of higher order schemes and turbulence modeling makes 
reasonably accurate simulations of extremely complex flowfields now possible. It is important to 
point out that if CFD codes are to be incorporated into a design process, the need for 
experimental validation is a necessity, a point that is sometimes under-emphasized. 4 

One of the advantages of CFD simulations over standard experimental data sets is that all 
flowfield properties are available everywhere in the flowfield so that any desired quantity can be 
extracted from the CFD results. Both detailed information about particular flowfield features and 
global information, such as conservation consistency checks and overall performance quantities, 
can be calculated from the numerical data set. 4 

The design of fuel injection configurations like the ones considered in the scramjet 
engines rely on the understanding of extremely complex flowfields with features such as: 
turbulent mixing, combustion reactions, strong three-dimensionality, compressibility effects, 
boundary layers and flow separations. The ability to study such flowfields experimentally has 
been limited because of the difficulty in gaining access into the harsh, yet easily perturbed 
environment of a supersonic flowfield. Understanding the geometrical complexity of such flows 
has been limited by the fact that most traditional quantitative measurement techniques are 
pointwise techniques. Theoretical studies in the past have been stalled by the difficulties in 
dealing with the highly non-linear set of equations that describe the physics of these flowfields. 
Numerical approaches give reasonable, and often very accurate, solutions to problems such as 
these that are not solvable analytically. Recently, due to the development nonintrusive laser 
diagnostics measurement techniques and to advances in the field of computational fluid 
dynamics, the ability to study such flowfields has been greatly improved. 4 
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Recent advances in numerical simulation techniques and theoretical studies offer new 
descriptions of the physical phenomena in high-speed flows including mixing and combustion 
interactions and the onset and development of instabilities and voticity, both in non-reacting and 
chemically reacting flows. Although numerical simulations offer much physical insight and a 
great level of detail, the experimental validation of both the theoretical models and the numerical 
algorithms is largely lacking due to difficulties in measuring essential parameters in these 
flowfields, such as velocity, temperature and pressure characteristics. Furthermore, there are 
numerous modeling issues, which are largely unresolved and need insight from experimental 
observation. Existing experimental facilities are, in general, dedicated to study of only a limited 
range of topics because of extreme conditions of high enthalpy flow, including an overview of 
the demands of availability of existing facilities. Certain facilities, capable of reproducing high 
enthalpy conditions, can operate for only short time durations, insufficient to achieve stable 
thermal conditions. Other existing facilities are, in general, limited in terms of their experimental 
flexibility. 6 
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2. THEORETICAL MODEL AND COMPUTATIONAL PROCEDURE 

The present computational study incorporates two cases with close configuration aiming 
at studying the effect of mixing schemes on the combustion characteristics in supersonic flow 
fields. In the first case study gaseous hydrogen is used as the main and only fuel in the system. 
While in the other case study gaseous hydrogen is used as a pilot energy source for the 
combustion of liquid kerosene. Both cases are treated three-dimensionally and both of them do 
have a rearward-facing step down stream of the test section entrance. The geometry of the 
configurations used in this computational study is presented in Figs. 2. 1-2.6. The flow 
configurations are matched with those employed in the experiments of Segal et al. 7 and Owens 
et al. 39 respectively. Figure 2.1 depicts schematically the configuration of the first case, which 
shows that the flow in the vicinity of the fuel injection stations is characterized by the supersonic 
flow of Mach number 2.0 over a rearward-facing step with transverse sonic gaseous H 2 injection 
downstream of the step. The length of the test section is 21.3 cm and the cross section is 2.54 by 
3.8 cm. The step height (H) is 5 mm in the upper wall. Hydrogen is injected normally to the 
incoming air flow direction from two stations, 3 and 7 step heights downstream of the step. The 
experimental injection diameters are 1 mm each. For the sake of facilitating the computational 
procedures and saving CPU time in the present computational study, the circular injection 
diameters were simplified to be square ones having the same cross section area. Because of the 
symmetry in the configuration under study, only one half of the physical domain was simulated 
with the plan of symmetry located at the centerline of the injection slots. Figure 2.2 presents a 3- 
D scheme for this configuration with the location of the plan of symmetry. Table 2.1 includes 
the flow conditions used in this case study. 
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In the second case study a vitiated airstream with 16% 0 2 mole fraction is directed to a 
test section entrance with a square cross section of 2.54 x 2.54 cm 2 . The test section, shown in 
Fig. 2.3, featured a 10-mm-high (H) rearward-facing step as the main flameholding mechanism. 
The step included a stagnation chamber from which the pilot flame is generated, parallel to the 
main flow via three, 1-mm-diameter orifices equally spaced in a transverse direction. Liquid 
kerosene (C 12 H 26 ) is injected normal to the incoming supersonic air flow through a hole of 
variable diameter of 0.5 and 1.0 mm. Removable, rearward-facing wedges were attached to the 
test section wall 2.8H downstream of the step with wedge length of 1.7H. The onset of thermal 
choking is delayed by diverging the test section starting 7H downstream of the step with a 3-deg 
half angle. Figure 2.4 presents the three injection configurations to be evaluated in this study. 
Configuration 1 featured a generic rectangular, rearward-facing step. Configuration 2 featured 
this step with either a 15 or 30-deg wedge installed 2.8H downstream of the step. Configuration 
3 featured the rearward-facing step with beveled-edge step which is shown in Fig. 2.4-c. The 
beveled step has the same height (10 mm) as the rectangular step, but the edges parallel to the 
main flow are beveled to facilitate vortex-enhanced mixing. 

Figures 2.5a-2.5c show schematic of the combustion flowfield for the three 
configurations. The air flow underwent expansion around the rearward-facing step where 
gaseous hydrogen is injected parallel to the main flow. When the pressure rise caused by 
hydrogen combustion was large enough, the expansion wave around the step turned into a 
compression shock. At five step heights downstream of the step, liquid kerosene is injected 
normal to the main flow creating a separation shock upstream of the jet. The jet then broke up 
downstream of the injection site where it vaporized and mixed with the main flow. Inclusion of 
the wedges generates the shock waves shown in Fig. 2.5-b, and the beveled step included 
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shedding of the counter rotating vortices as shown in Fig. 2.5-c. Each of these mechanisms 
interacts in a different way with the fuel jet. In all three-injection configurations, symmetry at the 
half width-plan exists. Therefore, only one half of the physical domain can be simulated. Figure 
2.6 shows the three gaseous-hydrogen-pilot injection holes at the base of the step with the 
symmetry plan located at the centerline of the middle hole. Table 2.2 includes the flow 
conditions used in this case study. 

The flow field to be analyzed is governed by the Reynolds averaged compressible form 
of the Navier-Stokes equations. These equations are closed using the two-equation k-e 
turbulence model. The equations are all-similar and can be expressed into the following common 
form: 

where the corresponding values for the exchange coefficient r<t> and the source or sink term So 
for any variable are summarized in Ref. [40], The Reynolds stresses putp and putp are related 
to mean flow quantities via K-e turbulence models. The transport equations for this model are 
available in Ref. [40]. 

The finite-volume method incorporated with the quadratic upwind interpolation scheme 
(QUICK) is employed to integrate and interpolate the partial differential equations represented 
by the common Eq. (2.1). The SIMPLEC, a variant of SIMPLE method developed by Patankar 
[41] was used to obtain the numerical solution for the discritized equations. The computer 
program used an Additive-Corrective Multi-Grid procedure to accelerate the convergence of the 
standard line-by-line solver. In addition to solving the basic transport equations on the original 
fine-grid level, multi-grid technique solve an equivalent equation set on each coarse grid level, 
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transferring correction onto the fine grid levels in an effort to achieve global balance and reduce 
the overall solution error. The use of multi-grid scheme can greatly reduce the number of 
iterations and CPU time required obtaining a converged solution. The convergence criterion was 
set so that the normalized residuals are less than 10' 6 at any grid node and for any <p equation. 
The full detailed description of the FULENT code flow solver is provided in Appendix A. 

For each case study the computations were performed on different grid sizes and the grid 
nodes were reasonably distributed and arranged to ensure that region of important influence on 
the flow field are adequately resolved. In general, the axial grid is clustered at the step and the 
fuel injection stations while the radial grid is concentrated within the jets, in the jets’ interaction 
zones and towards the walls. 

For the first case study a grid with a size of 209x34x23 is used. The aspect ratio was 
evaluated over the whole domain to make sure that the grid is fine enough to capture the flow 
features. 

Work in building the grid and setting up the simulation procedure for the second case 
study (pilot hydrogen with kerosene combustion case) is undercurrent investigation. 
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3. BOUNDARY AND INITIAL CONDITIONS 

No-slip boundary conditions are used along the combustor walls as well as along the step 
boundaries for all cases under study. The walls are assumed adiabatic requiring normal 
derivative of temperature to vanish. Along the supersonic inflow boundaries uniform conditions 
are used. Along the supersonic outflow boundaries, non-reflective boundary conditions are used 
where the boundary values are found by linear extrapolation from the interior points. 

Initial guessed values for the most important flow parameters such as temperature, 
velocity components, heat capacity (C p ), and mass fractions are assumed for the solution to start 
with. This process is a very important step in speeding the solution convergence. The inlet air 
values of these mentioned flow parameters are used in patching the whole flow field. The regions 
around the injection ports are patched using the corresponding values at the injection inlets. The 
patching process is explained in detailed in section 7.5.1 in Appendix A. 
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4. 


RESULTS AND DISCUSSION 


In this section some of the data obtained while running the first case (sonic gaseous 
hydrogen normal fuel injection) are presented. A global overview on the plan of work is 
presented after which the some figures are presented as part of the available results. 

This case was first simulated two-dimensionally and big differences between the 
experimental and numerically obtained wall pressure distribution were observed. Because of the 
three-dimensional nature of the normal injection it was necessary to guide the work towards the 
three-dimensional simulation. 

Following the step-by-step technique in the simulation approach, the cold flow field is 
solved first getting a converged solution. Thereafter, this cold flow converged solution is used as 
an initial solution for the simulation of the combustion flow field. A rough grid of 161x34x23 
was used in the beginning and converged solution was satisfied. 

The flow field was simulated first without normal fuel injection for just solving the 
supersonic flow over the rearward-facing step. This converged solution was used to calculate the 
flow field with normal fuel injection. A converged solution was obtained for such a case. 

The effect of grid refinement was tested by using more refined grid around the injection 
stations and in the far domain. A bigger grid with size of 209x34x23 was used and the cold flow 
field with fuel injection was simulated. The “Power Law Scheme” was used to descritize the 
governing equations. 

The effect of the descritization scheme was tested by switching on one of the higher order 
schemes (both the 2 nd order and quadratic schemes were used). The higher order schemes 
required more time to achieve the convergence criterion. Some illustrative results obtained are 
presented in this section. 
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Figure 4. 1 presents the normalized upper-wall pressure distribution along the combustor 
length for both “no-injection” and “with injection” cases. The averaged pressure at each axial 
location was normalized to the combustor entrance conditions. Both results were obtained using 
the smaller grid size of 161x34x23. It can be concluded that the pressure distribution with fuel 
flow and no combustion was essentially identical to those with no fuel flow. This conclusion is 
well matched with the experimental findings reported by Segal et al. 7 

The effect of grid refinement is shown in Fig. 4.2. The normalized upper-wall pressure 
distribution along the combustor length is plotted for both coarse and fine grids. A noticeable 
difference is clear especially just downstream of the step and in the far domain. This led us to 
discuss the effect of the descritization scheme. A higher order descritization scheme (2 nd order 
scheme) is compared to the power law scheme in Fig. 4.3. In this figure the normalized upper- 
wall static pressure distribution using both schemes are compared with the experimental values. 
The higher order scheme shows more comforting trend, especially just downstream of the step. 

The combustion flow field simulation is now underway using the higher order scheme 
solution as an initial solution. 
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5 . 


CONCLUSIONS 


Numerical investigation is conducted to study the effect of the different injection schemes 
on the mixing and thereupon the combustion of both gaseous hydrogen and liquid kerosene as 
two different alternative fueling systems in the SCRAM jet engine. The normal injection mode 
presents good potential towards increasing the mixing efficiency of the fuel and the supersonic 
incoming air. This effect helps in increasing the fuel residence time in the combustor therefore 
overcoming the difficulties of maintaining stable and efficient combustion in such very high 
speed flow fields. Injecting gaseous hydrogen parallel to the incoming air flow from the base of 
the step to form a pilot flame can help in maintaining the kerosene flame stability. It could be 
found that the maximum kerosene equivalence ratio that can maintain stable flame is mainly a 
function of the pilot hydrogen equivalence ratio. The flame stability limit can be increased by 
injecting more pilot hydrogen to some extent. The wedges increased the mixing efficiency by 
creating stream wise vortical structures. 

Work is still underway to complete simulating the different mixing schemes. 
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Table 1.1 Thermodynamic heats of reaction 


Hydrocarbon 

Reaction products 

Affi, 

Btu/lb 

Methylcyclohexane 

Toluene, Hi 

940 


Benzene, Hi, CHt 

710 


Cyclohexene, CH* 

286 


Cracked products 

-60* 

n-Oct&ne 

Octene, Hi 

491 


Octadiene, Hi 

ca 900 


StyTene, Hi 

1425 


Ethylbenzene (EB), Hi 

957 


Xylene, Hi 

915 


Tolueoe, Hi, CH 4 

718 


Benzene, Hi, C 1 H 1 

804 


Cracked products 

258 

ft-Dodecane 

Dodecene 

330 


Dodecadiene, Hi 

595 


Aromatics, Hi 

ca 640 


Cracked products 

- 110 - 

n-Hexadec? n© (cetane) 

Cracked products 

131 

Dicyclohexyl 

Benzene, Hi 

1038 


Toluene, pentane 

325 


Cyclohexane (CH) 

-94* 


Hexane 

-322* 


MCH, pentane 

-231* 


Light gas 

22 


Diphenyl, Hi 

1080 


Phenylcyclohexane, Hi 

540 

Decaiin 

Naphthalene, Hi 

950 


Tetralin, Hi 

670 


CH + olefin 

143 


MCH + olefin 

143 


DMCH + olefin 

143 


DECH 

- 110 * 


Benzene + olefin 

800 


Toluene + olefin 

800 


EB *4* olefin 

800 


DEB 

-700 


Alkyl aromatics 

-400 

Ethylbenzene 

Styrene 

509 


Toluene, CIL 

-273 


Benzene, CiH« 

-180 


Benzene, CiH* 

405 

JP-7 (F-71), 1200°F 

60% cracked, 50% gas 

310* 


• Exothermic. 

* Experiment*!. 
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Table 2.1 Free stream and injection nominal operating conditions for the transverse sonic 
fuel in jection case 


Ps)air (atm.) 

0.5 

Ttotaj)air (K) 

741 

Mach No.)air 

2.0 


1.4 

M)«r 

28.9 

Tstattc)fud (K) 

1500 

Psuiic)fuci (atm.) 

0.5 

Mach No.)fiiel 

1.0 


1.4 


i — 

K is the gas constant. 

*• 

M is the molecular weight 


Table 2.2 Fee stream and injection nominal operating conditions for the piloted-energy 
kerosene combustion case 


Mach No.)air 

1.8 

P totaI) air (Pa) 

430 

Ttotal)air (K) 

1000 

Equivalence Ratio ) h 2 

0.02 

Equivalence Ratio) k 

0.325 

Ttotal)k (K) 

298 and 420 
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vehicle underbody 



Fig. 1.1 Schematic of a generic scramjet engine configuration. 
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Fig. 1.2 Airframe-Integrated supersonic combustion ramjet (scramjet). 


MIXING CONTROLLED COMBUSTOR CONCEPT 



Fig. 1.3 Mixing controlled combustion concept. 
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UNDER-EXPANDED JET INJECTED NORMALLY 
INTO A SUPERSONIC FREE STREAM 
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Fig. 1.8 Typical flowfield features of under-expanded jet injected normally 
into a supersonic free stream. 




Unswept Ramp 



Swept Ramp 


Fig. 1.9 Ramp injector geometries. 


side view 
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ramp-generated varticie* 


Fig. 1.10 Ramp injector flowfield features 
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Fig. 1.12 Typical deposition behavior in air-standard jet fuels and model compounds. 
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Fig. 1.13 The energy per unit mass for various fuels. 








Fig. 1.15 The energy per unit volume for various fuels 





Fig. 1.16 Schematic diagram showing the complex delivry, distribution, and storage processes for 

cryogenic fuels. 
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Fig. 1.17 Losses of cryogenic liquid fuels dring refill and tank to tank transfer processes. 



Fig. 1.18 Heat-sink requirement as a function of flight speed. 
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Fig. 1.19 Possible modes of application of cooling for hydrocarbon fuels. 
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Fig. 1.21 Carrier-based hypersonic aircraft, inboard profile (to scale). 
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Fig. 2.1 Configuration of the sonic gaseous hydrogen injection case. 
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Fig. 2.2 3-D Configuration for the sonic gaseous 
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Fig. 2.3 Test section schematic for the kerosene 
combustion case. 


Fig. 2.4 Mixing Schemes used in the piloted-energy 
kerosene combustion case: a) baseline, 
b) shock-induced wedges, and c) beveled step 
for vortex-induced mixing. 
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Fig. 2.5 Combustion flowfield schematic: a) baseline, 
b) baseline with wedge, and c) beveled-step configurations. 
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Fig. 2.6 The pilot-hydrogen injection holes built in the step 
base of the piloted-energy kerosene combustion case. 


3D-Two-Hole Vertical Sonic H2 Injection 

Cold Flow with Total Air Temp.= 741 K 
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Fig. 4.1 Effect of injection on the flowfield of sonic gaseous hydrogen injection case. 



3D-Two-Hole Vertical Sonic H2 Injection 

Cold Flow with Total Air Temp.= 741 K 
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Normalized Axial Distance (X/H) 

lig. 4.2 Effect of grid refinement on the flowfleld of sonic gaseous hydrogen injection 





3D-Two-Hole Vertical Sonic H2 Injection 
Cold Flow with Total Air Temp.=741 K 


00 



(sj/d) 3Jnssajd paznBuuo^ 


Fig. 4.3 Effect of descritization scheme on the flowfield of sonic gaseous hydrogen injection case. 
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(Attached pages 1-67) 
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CFD “FLUENT” CODE 


1. INTRODUCTION 

FLUENT is a general purpose computer program for modeling fluid flow, heat transfer 
and chemical reaction. Using FLUENT you can quickly analyze complex flow problems even 
1 you do not have prior expertise in computational fluid dynamics or computer programming. 

FLUENT enables you to apply computer simulation methods to analyze and solve your 
practical design problems. 

FLUENT incorporates up-to-date modeling techniques and a wide range of physical 
models for simulating numerous types of fluid flow problems. These are accessible to you 
through an interactive graphical user interface for problem definition, computation and 
graphical postprocessing. When required, FLUENT can also be customized to your specific 
modeling needs and/or interfaced to your in-house CAD system. 

Other CFD Programs from Fluent Inc.: In addition to FLUENT, Fluent Inc provides 
other CFD programs: NEKTON, FLUENT/UNS, and RAMPANT. NEKTON is a finite 
element based program that uses the spectral-element method to solve fluid flow and heat 
transfer in complex and/or deforming geometries including non-Newtonian and free surface 
phenomena. NEKTON is ideally suited for prediction of materials processing and coating 
flows. FLUENT/UNS and RAMPANT are state-of-the-art computer programs for modeling 
fluid flow and heat transfer in complex geometries. FLUENT/UNS and RAMPANT provide 
complete mesh flexibility, solving your flow problems with unstructured meshes that can be 
generated about complex geometries with relative ease. Supported mesh types include 2D 
triangular/quadrilateral, 3D tetrahedral/hexahedral/pyramid/wedge, and mixed (hybrid) 
meshes. FLUENT/UNS and RAMPANT also allow you to refine or coarsen your grid based on 
the flow solution. RAMPANT utilizes a finite volume method tuned for highly compressible 
flows and FLUENT/ UNS utilizes a pressure-based finite volume method for incompressible 
and mildly compressible flows. FLUENT and FLUENT/UNS are targeted at similar 
applications, they differ only in mesh type (FLUENT accepts structured quad/hex meshes 
whereas FLUENT/UNS accepts meshes of all types, structured or unstructured) and range of 
physical models available. 

1.1 Basic Program Capabilities 

FLUENT can model a wide range of physical phenomena, including: 

• 2D/3D geometries in Cartesian, cylindrical or general curvilinear coordinates 

• steady state or transient flow 

• incompressible or compressible flow 

• laminar or turbulent flow, using a suite of turbulence models provided 

• coupled conduction/convection heat transfer (including both free and forced 
convection) 

• radiation heat transfer 
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• Mixing of chemical species, including optional multicomponent diffusion models 

• reaction of chemical species using either finite rate, equilibrium, or "mixed is 
burned" chemistry models 

• surface reactions, including chemical vapor deposition (CVD) 

• NO* formation and soot formation in combustion systems 

• multiphase flows in which two or more continuous phases are present 

• free surface flows with complex surface shapes 

• solidification/melting of liquids, with free surface effects and continuous casting 
models 

• deforming geometries using moving mesh models (e.g., valve closing) 

• sliding/rotating geometries using sliding mesh models (e.g., rotor/stator 
interactions) 

• temperature and composition dependent fluid/material properties 

• laminar flow of non-Newtonian fluids 

• flow through porous media, including thermal effects of the solid media 

• one-dimensional fan/heat-exchanger performance models 

• dispersed second phase of particles/bubbles/droplets, including 

- Lagrangian trajectory calculations with stochastic tracking to account for the 
effects of turbulence, 

- inert heating or cooling of the dispersed phase, 

- evaporation of the liquid droplets into the continuous phase, 

- coal particle devolatilization and char burnout, 

- heterogeneous surface reactions 

- coupled mass, momentum, and heat transfer between the dispersed and 
continuous phases 

The size and scope of your problem (in terms of the number of computational nodes 
and the number of chemical species and/or reactions) is limited only by your available 
computer memory. 

Applications: The physical models listed above enable you to apply FLUENT to a 
wide range of applications, such as: 

• chemical and process engineering component design 

• combustion design and engineering, including gaseous combustion, liquid fuel 
combustion, and coal combustion 

• aerodynamic design 

• electronics cooling, manufacture, and design 

• power generation 

• heat transfer operations 

• materials processing 

• chemical vapor deposition (CVD) 

• spray drying or cooling 
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• gas cleaning or particle classification 

• architectural design (internal and external air flow) 

• fire research (open fires, fires in buildings) 

• manifold design 

• particulate deposition or fouling 

• pollution control 

• turbomachinery component design 

Solution Technique: FLUENT models this wide range of phenomena by solving the 
conservation equations for mass, momentum, energy, and chemical species using a control 
volume based finite difference method. The governing equations are discretized on a 
curvilinear grid to enable computations in complex/irregular geometries. A nonstaggered 
system is used for storage of discrete velocities and pressures. Interpolation is accomplished 
via a first-order, Power-Law scheme or optionally via higher order upwind schemes. The 
equations are solved using SIMPLE-like algorithms with an iterative line-by-line matrix solver 
and multigrid acceleration. 

Boundary Conditions: You can define the unique conditions that describe your 
problem via a wide variety of boundary conditions, including: 

• flow inlets defined via velocity 

• flow inlets and/or exits defined via pressure 

• flow exits defined via extrapolated/fully-developed outlet conditions 

• stationary or moving/rotating walls 

• slip (frictionless) or non-slip walls 

• isothermal walls or walls with a specified heat flux or external heat transfer 
condition 

• conducting solids 

• periodic (fully-developed flow) boundaries 

• symmetry boundaries 

• cyclic boundaries 

1.2 Program Structure 

FLUENT is part of the Fluent Suite, which includes the following products: 

• FLUENT, RAMPANT, FLUENT/UNS, and NEKTON, the solvers. 

• prePDF, the preprocessor for modeling PDF combustion. 

• GeoMesh and preBFC, the preprocessors for geometry modeling, block-structured 
mesh generation, unstructured triangularsurface mesh generation, and 2D 
unstructured triangular mesh generation. 

• TGrid, the 2D triangular and 3D tetrahedral mesh generator. 

grid filters for CAD/CAE packages such as I-DEAS surface and volume grid 
generators. 
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Figure A.l shows the organizational structure of these components. (TGnd does not 
appear in the diagram because it does not create grids for use in FLUENT.) 

preBFC and GeoMesh are used to define the geometry and create a structured grid. 
preBFC is intended for relatively simple 3D geometries, whereas GeoMesh has considerably 
more geometry modeling and mesh generation capabilities. 

Quadrilateral and hexahedral grids created with GeoMesh can be saved in FLUENT 
format and read directly into FLUENT. Quadrilateral and hexahedral grids created with 
preBFC are automatically written in FLUENT format and can be read directly into FLUENT. 
Grid generation utilities for simple Cartesian and cylindrical-polar meshes are also available 
within FLUENT. 

In addition, volume grids generated in ICEMCFD, I-DEAS, PATRAN, ANSYS, and 
GREDGEN can be imported into FLUENT. 

Once a grid has been read into FLUENT, all remaining operations are performed within 
the solver. These include setting boundary conditions, defining fluid properties, executing the 
solution, refining the grid, and viewing and postprocessing the results. 

2. GETTING STARTED WITH FLUENT 

FLUENT features a graphical user interface (GUI). Most operations, including file 
reading/writing, problem setup, calculation, and postprocessing, can be accomplished via the 
GUI. The menu-driven text interface is also available for these and additional operations. In 
this chapter, the use of the GUI and the text interface is demonstrated by setting up and solving 
a sample problem. 

2.1 Overview of Using FLUENT 

FLUENT and Preprocessors: FLUENT is designed to run as a stand-alone program, 
or in conjunction with GeoMesh or preBFC. Geometries that can be modeled with Cartesian 
or cylindrical coordinates in either 2D or 3D may be set up and solved using FLUENT alone. 
For other analyses in which body-fitted coordinates are required, GeoMesh or preBFC can be 
used to set up the geometry and grid. 

Planning Your CFD Analysis: When you are planning to solve a problem using 
FLUENT, you should first give consideration to the following issues: 

• Definition of the Modeling Goals: 

What specific results are required from the CFD model and how will they be used? 

What degree of accuracy is required from the model? 

• Choice of the Computational Model: 

How will you isolate a piece of the complete physical system to be modeled? 

Where will the computational domain begin and end? What boundary conditions 
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will be used at the boundaries of the model? Can the problem be modeled in two- 
dimensions or is a three-dimensional model required? 

• Design of the Grid: 

Will your model use a Cartesian grid or a Body-Fitted Coordinate (BFC) grid? 
What grid topology will work best in Body-Fitted Coordinates? What degree of 
accuracy will be required in each region of the domain? Where are the flow details 
that will require a refined grid? 

• Choice of Physical Models: 

Is the flow laminar or turbulent? Which turbulence model should be applied? Is the 
flow steady or unsteady? Is heat transfer important? Will buoyancy effects be 
significant? Will you treat the fluid as compressible or incompressible? Are there 
other physical models that should be applied? 

• Determination of the Solution Procedure: 

Can the problem be solved simply using the default solution parameters in 
FLUENT? Can convergence be accelerated via a more judicious solution 
procedure? Perhaps the problem involves "highly coupled" physics, in which the 
governing equations have a high degree of interdependency. Will you need a step- 
wise solution procedure to obtain a stable solution? How long will the problem take 
to converge on your computer? 

Careful consideration of these issues before beginning your CFD analysis will 
contribute significantly to the success of your modeling effort. 

Problem-Solving Steps: Once you have determined the important features of the 
problem you want to solve, the basic procedural steps you will follow are those shown below. 
Steps that are optional are noted as such. 

1. Create or import the model geometry and grid 

2. Choose the basic equations to be solved (e.g., enthalpy, species, turbulence transport) 

3. Identify additional models needed (fans, porous media, special boundary conditions, 
species transport or chemical reaction, etc.) 

4 . Specify the boundary conditions 

5. Specify the fluid properties 

6. Set up a dispersed phase (optional) 

7. Adjust the solution control parameters (optional) 

8. Calculate a solution (fluid phase and/or dispersed phase) 

9. Examine the results 

10 . Save the results 

11. Consider revisions to the numerical or physical model 
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3. 


GEOMETRY SETUP/GRID GENERATION 


3.1 Overview of Geometry and Grid Definition Options 

FLUENT employs what could be termed a gnd-based geometry, in which the geometry 
of your model is determined by control volumes defined by the grid. You can choose to use 
either of the following: 

• Curvilinear body-fitted coordinates (BFC) in which the grid lines are determined by 
a coordinate system that conforms to the geometric boundary of your model. 

• Cartesian or cylindrical coordinates in which the grid lines are aligned with the 
cartesian (x, y, z) or cylindrical (r, 0, z) coordinates and may not conform to oblique 
or curved boundaries of your physical geometry. 

In addition, you can choose one of several techniques to create the grid, using preBFC, 
GeoMesh, or FLUENT, or importing the grid from an external CAD package. These grid 
generation choices are discussed in more detail in this section. 

Body-Fitted Coordinates vs. Cartesian or Cylindrical Coordinates: Figure A.2 
depicts a problem that has been defined first with body-fitted coordinates and then with 
Cartesian coordinates. The approximation of the geometry via Cartesian or cylindrical 
coordinates is evident in the second example. Generally, you will choose to use body-fitted 
coordinates if your physical geometry is not approximated well by a simple Cartesian or 
cylindrical coordinate system and you are concerned with an accurate description of the flow 
and/or heat and mass transfer at the physical boundaries. Sometimes, however, you may 
choose the simplicity of a Cartesian or cylindrical coordinate system, even if the physical 
boundaries of the problem must be approximated by a "stair-step" pattern in the grid. Figure 
A.3 illustrates one such "stair-stepped" model in which the accuracy of the predicted flow field 
is unlikely to be severely impacted by the approximation of the curved boundary. This is often 
true when the flow field is dominated by internal shear layers and jet mixing as in the example 
shown. However, when the wall boundary layers are of great importance to the overall flow 
pattern, an accurate body-fitted description of the boundary is important, as illustrated in the 
diverging duct example shown in Fig. A.4. 

Preprocessor or FLUENT? If you choose to use body-fitted coordinates (BFC), you 
must define your problem geometry and the computational grid using preBFC or GeoMesh. If 
you choose to use Cartesian or cylindrical coordinates, you can define the problem geometry 
and grid in preBFC or you can perform the grid generation directly in FLUENT. Figure A. 5 
illustrates these alternate setup options. 

3.2 Geometry/Grid Requirements and Considerations 
3.2.1 Geometry/Grid Requirements in FLUENT 
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FLUENT imposes several geometry setup and grid construction requirements that you 
should be aware of before beginning your problem setup: 

• Axisymmetnc geometries must be defined such that the axis of rotation is the X- 
axis of the Cartesian coordinates used to define the geometry (Fig. A. 6). 

• Three-dimensional geometries that are rotationally cyclic must be defined such that 
the axis of rotation is the Z-axis of the Cartesian coordinates used to define the 
geometry (Fig. A. 7) when the Cartesian velocity formulation is used. The axis of 
rotation should be defined as the X-axis when the cylindrical velocity formulation is 
used. 

• Physical boundaries that are to be modeled as cyclic boundaries must be located at 
the grid planes I = 1 and I = NI where NI is the maximum I index in your grid. (See 
Figs. A. 7 and A. 8). Cyclic boundaries cannot be located on planes of constant J or 
K. For rotationally cyclic geometries, J should be the radial direction and K should 
be the axial direction. 

• Interior regions (e.g., regions of dead or wall cells that are surrounded by live cells 
on more than one edge) should, ideally, be at least two control volumes wide and 
must be at least one control volume wide in each coordinate direction. This 
restriction exists due to the logic used in FLUENT to store quantities associated 
with the wall surface. An interior wall region that is only one cell thick may be 
used but will not provide sufficient storage for correct reporting of surface 
quantities (pressures, heat fluxes, etc.) associated with the two live cell interfaces on 
either side of it. Furthermore, there are two cases in which an interior wall region 
that is only one cell thick should not be used because it can adversely affect the 
calculation: 

- The combination of one-cell-thick walls and a cut link for any variable should 
not be used. 

- One-cell-thick walls cannot be used with the VOF model. 

Care should be taken to observe these constraints when defining your problem 
geometry and computational grid. 

3.2.2 Choosing a Good Grid 

Accuracy and Stability: The grid is a discrete representation of the continuous field 
phenomena that you are modeling and the accuracy and numerical stability of your FLUENT 
simulations depend on your choice of grid. In other words, the density and distribution of the 
grid lines determines the accuracy with which your FLUENT model represents the actual 
physical phenomena. 
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Differencing Schemes: In FLUENT, the control volume method, sometimes referred 
to as the finite volume method, is used to discretize the transport equations. In the discrete 
form of the equations, values of the dependent variables appear at control volume boundary 
locations (cell faces). These values have to be expressed in terms of the values at the nodes of 
neighboring cells in order to obtain algebraic equations. This task is accomplished via an 
interpolation practice, also called a differencing scheme ". The choice of differencing scheme 
not only affects the accuracy of the solution but also the stability of the numerical method. 

The use of an interpolation assumption introduces a discretization error, which 
decreases as the grid spacing is reduced. In practice, it is necessary to obtain sufficiently 
accurate solutions by employing grids that are not excessively fine. One measure of the 
discretization error of a scheme is the so-called truncation error of a Taylor series 
approximation to the derivatives in the governing differential equations. Truncation error 
occurs due to the approximate nature of the finite difference representation and can be reduced 
by increasing the order of the interpolation method. You should note, however, that for a given 
grid, higher order schemes (i.e., those with smaller truncation errors) do not necessarily lead to 
smaller solution errors. The discretization error is also related to how well the interpolation 
assumption represents the real physical phenomena. For example, for a one-dimensional 
situation, in which diffusion plays a dominant role in the conservation equation, a linear profile 
for the dependent variable is expected. However, if both convection and diffusion are present, 
but convection plays a dominant role, a linear profile introduces large errors since an 
exponential profile is expected. 

Numerical Diffusion: In addition, a dominant source of errors in multi-dimensional 
situations is the so-called false diffusion or numerical diffusion. The term false diffusion is 
used because its effect on a flow field is analogous to that of increasing the real diffusion 
coefficient. False diffusion is only noticeable when the real diffusion is small, that is, when the 
situation is convection-dominated. The phenomenon of false diffusion arises due to one- 
dimensional interpolation practices being employed in a multi-dimensional situation. 
Therefore, this source of error occurs when the flow is oblique to the grid lines and when there 
is a nonzero gradient of the dependent variable in the direction normal to the flow direction. 

The ideal differencing scheme, therefore, has the following attributes: 

• low truncation error 

• low numerical diffusion 

• physically meaningful variation 

The first and second attributes directly affect the accuracy of the solution. The accuracy 
can be increased by increasing the number of grid lines. 

Differencing Schemes in FLUENT: The default differencing scheme used in 
FLUENT is the so-called power-law scheme. This scheme is derived from the exact analytical 
solution to the one-dimensional convection-diffusion equation. The power-law scheme is very 
stable and gives physically meaningful (bounded) solutions but, in certain situations, is 
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susceptible to numerical diffusion effects. These effects, as mentioned earlier, are maximum 
when the flow is aligned at 45 degrees to the grid lines and there are significant gradients in the 
direction normal to the flow. Higher-order methods which are less susceptible to numerical 
diffusion, but also less stable compared to the power-law scheme, are also available in 
FLUENT. 

Solutions in which truncation error and/or numerical diffusion are significant are 
termed grid-dependent. If the grid is refined until the solution no longer varies with additional 
grid refinement, you have achieved a gnd-independent solution. When you are interested in 
quantitative accuracy, you should ensure grid independence of the solution in this way. 

Resolution of Flow Details: The overall number of grid lines to be used in a FLUENT 
simulation is usually determined by noting the locations and sizes of expected flow features 
such as shear layers, separated regions, boundary layers, and mixing zones, and then choosing 
a sufficient number of grid lines to resolve these important flow details. Typically, inlets and 
outlets should be defined by a minimum of two or three cells and are best represented via 10 to 
12 cells. Where large gradients are expected, as in shear layers or mixing zones, the grid should 
be fine enough to minimize the change in the flow variables from node to node. Separated 
regions should include at least 5 or 6 cells across the separation and are best represented with 
more cells. In general, no flow passage should be represented by fewer than 3 or 4 cells. 
Figure A. 9 illustrates the resolution of the flow using two different grids in a simple driven 
cavity. 


Grid Spacing Near Walls: The spacing between a wall and the adjacent grid line can 
impact the accuracy of the computed shear stress and heat transfer coefficient at the wall. This 
is particularly true in laminar flow where the grid adjacent to the wall should obey 





(A.l) 


Where y p 

Uoo 

v 

x 


distance to the wall from the adjacent node 
free-stream velocity 
kinematic viscosity of the fluid 

distance along the wall from the starting point of the boundary layer. 


Equation (A.l) is based upon the computed boundary layer profile in laminar flow over 
a flat plate at zero incidence. This criterion ensures that the gradients at the wall (e.g., 5u/8n) 
are adequately approximated by the difference expression (e.g., (Up - U wa ii)/An). The situation 
is depicted visually in Fig. A. 10. 


Similar guidelines for the choice of the near-wall grid spacing in laminar flows can be 
derived based on the analytical solutions for fully-developed flows. For fully-developed 
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laminar flow in a pipe of radius R, you should choose — <0.1 and for fully-developed 

R 

laminar flow between parallel plates separated by a height H, you should choose — <0 05 

H 

These guidelines yield prediction of x w , within 10% of the theoretical values for fully 
developed flows. Placement of the grid lines closer than these guidelines improves the 
accuracy of your predictions. In developing flows, where the boundary layers are thinner, the 
near-wall grid spacing should be smaller than that based on these guidelines. 

Non-Uniform Grid Spacing: Often the guidelines noted above may give rise to a very 
large number of computational cells. One way to minimize the number of cells while 
maintaining a sufficient degree of accuracy in the solution is to use a non-uniform grid. In a 
non-uniform grid, the grid spacing is reduced in regions where high gradients are expected and 
increased in regions where the flow is relatively uniform. 

Certain guidelines should be followed in the generation of non-uniform grids. Most 
importantly, the rate of change of grid spacing should be minimized. Normally the spacing 
between adjacent grid lines should not change by more than 20% or 30% from one grid line to 
the next. This implies expansion factors between 0.7 and 1.3. This is an accuracy 
consideration, primarily impacting the accuracy of the diffusion terms in the governing 
transport equations. 

Cell Aspect Ratios: The aspect ratio of the computational cells is an additional issue 
that arises during the setup of the computational grid. While large aspect ratios may be 
acceptable in some problems, a general rule of thumb might be to avoid aspect ratios in excess 
of 5:1. This limit can be acceptably exceeded when the gradients in one direction are very 
small relative to those in the second direction, such as occurs in a fully developed pipe flow. 
Conversely, excessive aspect ratios can lead to stability problems, convergence difficulties, 
and/or the propagation of numerical errors. 

Grid Skewness: When you use Cartesian coordinates to define your model, the 
coordinate system guarantees that the grid lines intersect at 90 degree angles. When you use 
body-fitted coordinates, the grid lines you generate may not be orthogonal. While some degree 
of non-orthogonality is allowable, and is accounted for in the solution process, your 
computational grid should maintain grid intersection angles close to 90 degrees whenever 
possible. Figure A. 11 depicts two examples of grids with excessive skewness (departure from 
orthogonality). Whenever possible you should limit the skew in your grid to less than 45 
degrees. 


3.3 Grid Lines, Nodes, Control Volumes, and Cell Types 

In order to understand the relationship between the computational grid and geometry, 
you need to understand the grid definition and cell numbering system used by FLUENT. 
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Figure A. 12 illustrates this nomenclature which is described below. Note that, internally, 
FLUENT uses a different numbering system. See the user-defined subroutine chapter for 
details. 


Grid Lines and Nodes: Figure A. 12 shows solid lines which are grid lines (lines of 
constant I or lines of constant J). The grid lines intersect at node points, with the intersection of 
grid lines I and J forming node ( 1 ,J). 

Control Volumes: The grid lines define the boundaries of control volumes or cells. 
These control volumes are used when integrating the governing equations and form the basis of 
your computational model. Each control volume in your FLUENT model has a cell type. The 
cell type tells FLUENT whether the control volume is filled with fluid, or if the control volume 
defines a wall, inlet, exit, etc. 

Cell Center/Storage Locations: The cell center ( 1 ,J) is located at the geometric center 
of the control volume or cell (I,J). This cell center is the storage location for all dependent 
variables in FLUENT. 

Grid and Cell Numbering: Figure A. 12 shows the control volume (1,J) which is 
located between grid lines 1-1 and I and grid lines J-l and J. That is, the grid point (I,J) is 
located at the upper right comer of control volume or cell (I,J). In 3D, grid point (I,J,K) is 
located at the upper right, front comer of cell (1,J,K). 

The relationship between grid numbers and cell numbers is clarified in Figure A. 13. 
Figure A. 13a shows a matrix of 6 grid lines in each computational direction, i.e., a 6 x 6 grid. 
The grid points (1,1) and (6,6) are indicated in the figure. These grid lines define a block of 5 
control volumes in each coordinate direction. The cell numbering is. noted in Figure A. 13b, 
with the first cell in the lower left comer denoted as cell (2,2) and the last cell in the upper right 
comer denoted as cell (6,6). 

Fictitious Boundary Cells: FLUENT adds an additional ring of cells or control 
volumes around the perimeter of your computational model. These control volumes are used to 
define boundary conditions at the edge of the physical domain. Thus, the complete FLUENT 
description of the example shown in Fig. A. 13 becomes that of Fig. A. 14, in which there is a 
matrix of 7 x 7 control volumes. The first cell in the lower left comer is cell (1,1) and the final 
cell in the upper right comer is cell (7,7). Both of these cells are part of the group of fictitious 
boundary cells surrounding the matrix of 5 x 5 physical control volumes. 

4. BASIC PHYSICAL MODELS FOR FLOW AND HEAT TRANSFER 

This chapter describes the basic physical models that FLUENT provides and the input 
commands that you use to define and use them. 

4.1 Selecting Physical Models and Defining the Problem Scope 
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When you are beginning the setup of a new FLUENT model you must begin by 
defining the problem scope in terms of the basic physical models that should be applied 
FLUENT’S default physical models consist of laminar, isothermal, single-component 

flow. However, you can alter this problem scope, as described in this section, to include 
the following: 

• Cylindrical velocity formulation for 3D problems 

• Turbulent flow 

• Heat transfer 

• Radiation heat transfer 

• Periodic flow and heat transfer 

• Porous media 

• Compressible flow 

• Time dependence 

• Lumped parameter fan models 

• Lumped parameter heat exchanger models 

• Chemical species transport and/or reaction 

• Phase change 

• Multiphase flow 

• Swirling and rotating flows 

• Rotating reference frames 

• Sliding mesh 

• Deforming mesh 

• Non-Newtonian viscosity 

4.2 Turbulence Models 

In turbulent flows, the velocity at a point is considered as a sum of the mean 
(ensemble-averaged) and fluctuating components: 


u ‘ ~ u t +u , (A. 2) 

Substituting expressions of this form into the instantaneous momentum equations 
(and dropping the overbar on the mean velocity, u) yields the ensemble-averaged 
momentum equations: 
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Equation (A. 3) has the same form as the fundamental momentum balance with 
velocities now representing ensemble-averaged (or mean-flow) values. The effect of 
turbulence is incorporated through the "Reynolds stresses", pu,u,. FLUENT relates the 
Reynolds stresses to mean flow quantities via one of three turbulence models: 

• The k-E model 

• The RNG k-E model 

• The Reynolds Stress Model (RSM) 

4.2.1 Choosing the Appropriate Turbulence Model 

The process of selecting a turbulence model for a given turbulent flow problem is 
greatly facilitated when you have a good understanding of the salient features of the flow 
in question. Based on this understanding, you should then consider which model would 
be more suitable. To do so, you must know the capabilities and limitations of the 
individual models. This section provides you with an overview of the models and general 
guidelines that will help you choose the correct turbulence model for the flow you want 
to model. 

The standard k-E model proposed by Jones and Launder has been the workhorse 
of engineering turbulence models for more than two decades. It falls in the category of 
"two-equation" turbulence models based on an isotropic eddy-viscosity concept. As such, 
it is more universal than other low-order turbulence models such as algebraic ("zero- 
equation”) and one-equation models. Robustness, economy, and reasonable accuracy for 
a wide range of turbulent flows explain its popularity in industrial flow and heat transfer 
Simulations. It is a semi-empirical model, and the derivation of the model equations, 
including the various model constants, relies on phenomenological considerations and 
empiricism. 

The Renormalization Group (RNG) k-E model also belongs to the k-E family of 
models. The model equations in their RNG form are similar to those for the standard k-E 
model. There are major differences, however, between the RNG and standard k-E models. 
The RNG model was derived using a more rigorous statistical technique, and its model 
constants are derived "analytically". 

The Reynolds Stress Model (RSM) is the most elaborate turbulence model that 
FLUENT provides. Eschewing the isotropic eddy viscosity hypothesis, the RSM closes 
the Reynolds stresses by solving their transport equations (6 additional equations in 3D, 
in comparison with k-E models). As such, the RSM accounts for the history and transport 
of the Reynolds stresses in a rigorous manner. The effects of streamline curvature, swirl, 
and rotation are all directly accounted for by the transport equations for the Reynolds 
stresses. 

Major Differences Between the RNG k-E Model and the Standard k-E Model: 

The RNG k-E model differs from the standard k-E model in several important ways: 


13 



• While the standard k-E model is based on the traditional Reynolds-averaging 
technique, the RNG model is derived by a more rigorous statistical (scale- 
elimination) technique. 

• The RNG model has an additional term in its e equation which significantly improves 
the accuracy for rapidly strained flows. 

• The effect of swirl on turbulence is included in the RNG model, which enhances the 
accuracy for swirling flows. 

• The RNG theory provides an analytical formula for turbulent Prandtl numbers, while 
the standard k-E model uses user-specified, constant values. 

• While the standard k-e model is a high-Reynolds-number model, the RNG theory 
provides an analytically-derived differential formula for effective viscosity that 
accounts for low-Reynolds number effects. 

All the features listed above make the RNG k-E model more accurate and reliable 
for a wide class of flows than the standard k-E model. 

Examples of flows for which use of the RNG k-E model is particularly beneficial 
are as follows: 

• Separated flows and recirculating flows (e.g., backward- or forward-facing steps, 
sudden expansions, diffusers, bluff bodies, lift devices at a high incidence) 

• Flows in curved geometries and flows that are rapidly strained (e.g., curved ducts, 
bluff bodies, highly accelerating/decelerating flows, stagnation and detaching flows) 

• Time-dependent flows with large-scale organized structures (e.g., vortex shedding, 
shear-layer instability) 

• Heat transfer in low-Prandtl-number fluid flows (e.g., liquid metal flows) 

• Low-Reynolds-number or transitional flows (e.g., flows in an otherwise quiescent 
enclosure where the flow is turbulent in regions of limited extent, but is otherwise 
laminar) 

• Flows with streamwise vortices and secondary flows (e.g., secondary flows in curved 
ducts and transition ducts, horseshoe vortices around junctions, streamwise vortices 
behind aerodynamic/hydrodynamic bodies such as ground vehicles, under-water 
bodies, and airplanes) 

• Swirling flows (e.g., swirl combustors, cyclones) 

If you expect that the flow you are modeling exhibits any of the features listed 
above, it is recommended that you use the RNG k-E model. 

RNG k-E Model vs. the Differential Reynolds Stress Model: Although the 
RNG model provides substantial improvements over the standard k-E model as stated 
above, it should be understood that the RNG model is still based on the isotropic eddy- 
viscosity concept. Consequently, there is a limit to what the RNG model can offer when 
the anisotropy of turbulence has a dominant effect on the mean flow. Such cases include 
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highly swirling flows, flows with strong streamline curvature, stress-driven secondary 
flows, and evolution of streamwise vortices. Among these, for instance, stress driven 
secondary flows cannot be predicted by the RNG model. In the RSM, all of these effects 
are automatically accounted for in the transport equations for the Reynolds stresses. 
When your simulation requires rigorous modeling of these features, you should consider 
using the RSM. 


Turbulence Modeling in Swirling Flows: If your FLUENT model involves a 
highly swirling flow, you should use the Reynolds Stress Model or the RNG k-z model to 
obtain a more accurate flow prediction. Generally, these models provide better accuracy 
in modeling devices such as cyclone separators or swirl nozzles in which the Swirl 
Number, S, approaches or exceeds unity. The Swirl Number can be defined as the ratio of 
the axial flux of angular momentum to the axial flux of axial momentum: 


S = 


| rv e v.dA 
R^uv.dA 


(A.4) 


Conservation of angular momentum (r 2 Q = constant) creates a tendency in the 
fluid toward a free-vortex form as depicted in Fig. A. 15. When the angular momentum in 
the fluid is quite high, this form may not be captured by the k-e model due to the 
assumption of isotropy used in the calculation of the Reynolds stresses. If your FLUENT 
model includes highly swirling flow, you should consider applying the Reynolds Stress 
Model (RSM) or the RNG k-z model. 


Computational Effort: CPU Time and Solution Behavior: Due to the extra 
terms and functions in the governing equations and a greater degree of nonlinearity, 
computations with the RNG k-z model tend to take more CPU time than with the 
standard model. However, the difference is marginal (10-15% increase compared to the 
standard k-z model). 

Since diffusion has a stabilizing effect on the numerics, the k-z model, which 
usually gives larger turbulent diffusivity, is more stable. The RNG model is more likely 
to be susceptible to instability. However, you should understand that this is because the 
RNG model, where appropriate, yields a correct level of turbulent viscosity (which is 
usually smaller than that computed by the standard k-z model). The RNG model is more 
faithful in responding to any physical instabilities such as time-dependent turbulent 
vortex shedding. 

Compared with the k-z models, the RSM requires additional memory and CPU 
time due to the increased number of the transport equations for Reynolds stresses. 
However, efficient programming in FLUENT has reduced the CPU time per iteration 
significantly. On average, the RSM in FLUENT requires 50-60% more CPU time per 
iteration compared to the k-z models. It should be noted, however, that the inter-coupling 
among the Reynolds stresses and mean flow makes solution convergence relatively slow, 
and that 15-20% more memory is needed. 
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Enthalpy Sources Due to Reaction: Sources of enthalpy, S h , in Equation (A. 5) 
include the source of enthalpy due to chemical reaction: 
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(A. 9) 


where h. is the enthalpy of formation of species j’ and R is the volumetric rate of 
creation of species j’. 


4.4 Compressible Flows 


Compressibility effects are encountered in gas flows at high velocity and/or in 
which there are large pressure variations. When the gas flow velocity approaches or 
exceeds the speed of sound or when the pressure change in the system (Ap/p) is large, the 
variation of the gas density with pressure has a significant impact on the flow velocity, 
pressure, and temperature. Compressible flows create a unique set of flow physics for 
which you must be aware of the special input requirements and solution techniques 
described in this section. Figure A. 16 illustrates several examples of compressible flows 
computed using FLUENT. 

When to Use the Compressible Flow Model: Compressible flows can be 
characterized by the value of the Mach number 


M -u!c 


(A. 10) 


where c is the speed of sound in the gas: 




(A. 1 1) 


and v is the ratio of specific heats (Cp/C v ). 

When the Mach number is less than one, the flow is termed subsonic. At Mach 
numbers much less than one (e.g., M < 0.3 or so), compressibility effects are negligible 
and the variation of the gas density with pressure can safely be ignored in your flow 
modeling. When the Mach number approaches unity, however, compressibility effects 
become important. The flow will choke at Mach 1.0. When the Mach number exceeds 
1.0, the flow is termed supersonic. Supersonic flows may contain shocks and expansion 
fans which can impact the flow pattern significantly and which require compressibility in 
your FLUENT model. FLUENT provides a full range of compressible flow modeling 
capabilities for subsonic, transonic, and supersonic flows. 

Special Physics of Compressible Flows: Compressible flows are typically 
characterized by the total pressure P 0 (isentropic stagnation pressure) and total 
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temperature T 0 (isentropic stagnation temperature) of the now. These quantities can be 
related to the static pressure and temperature via the following relationships: 


T y — 1 

(A. 13) 

These unique relationships describe the variation of the static pressure and 
temperature in the How as the velocity (Mach number) changes under isentropic (loss- 
free, constant enthalpy) conditions. 

For a given pressure ratio from inlet to exit, for example. Equation (A. 12) can be 
used to estimate the exit Mach number (which would exist in a one-dimensional 
isentropic flow). For air, Equation (A. 12) predicts a choked flow (Mach number of 1 0) at 
an isentropic pressure ratio, P s /P 0 , of 0.5283. This choked flow condition will be 
established at the point of minimum flow area (e.g., in the throat of a nozzle). In the 
subsequent area expansion the flow may either accelerate to a supersonic flow in which 
the pressure will continue to drop, or the flow may return to subsonic flow conditions, 
decelerating with a pressure rise. When a supersonic flow is exposed to an imposed 
pressure increase, a shock will occur, with a sudden pressure rise and deceleration 
accomplished across the shock. 

Basic Equations for Compressible Flows: Compressible flows are described by 
the standard continuity and momentum equations solved by FLUENT, and you do not 
need to activate any special physical models (other than the compressible treatment of 
density as detailed below). The energy equation solved by FLUENT, Equation (A.5), 
correctly incorporates the coupling between the flow velocity and the static temperature, 
and should be activated whenever you are solving a compressible flow. In addition, you 
should activate the optional viscous dissipation terms in Equation (A.5), which become 
important in high-Mach-number flows. 

5. CHEMICAL SPECIES TRANSPORT AND REACTING FLOW 

FLUENT models chemical species transport and chemical reactions using the 
reacting flow models described in this chapter. You can model any number of chemical 
species and any number of chemical reactions, limited only by the memory available on 
your computer. This chapter describes the chemical species modeling options in 
FLUENT, the options that you use to define chemistry, and the computational techniques 
that can be employed to yield a successful simulation. 

Reacting Systems Handled by FLUENT: Chemical reactions in FLUENT may 
be of the following types: 
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• Fluid phase reactions that may involve multiple chemical reaction steps, equilibrium 
chemistry, gaseous combustion, NO, and other pollutant formation, etc. 

• Liquid fuel combustion in which fuel vapor is generated via evaporation of liquid 
droplets and a combustion reaction occurs in the gas phase. 

• Surface reactions (e.g. Chemical Vapor Deposition) in which the reaction occurs at a 
solid (wall) boundary. 

• Particulate reactions (e.g. coal combustion) in which reaction occurs by the evolution 
of a combustible gas and/or at the surface of a solid particle. 

This section describes how you can model these reacting systems. 


Reaction Modeling Options: FLUENT provides three reaction modeling approaches: 

1. Generalized Finite Rate Formulation 

This approach is based on the solution of species transport equations for reactants 
and product concentrations, with the chemical reaction mechanism defined by you. The 
reaction rates that appear as source terms in the species transport equations are computed 
from Arrhenius rate expressions or by using the eddy dissipation concept due to 
Magnussen and Hjertager. Models of this type are suitable for a wide range of 
applications including combustion simulation. 

2. Mixture Fraction/PDF Formulation 

In this approach individual species transport equations are not solved. Instead, 
transport equations for one or two conserved scalars (the mixture fractions) are solved 
and individual component concentrations are derived from the predicted mixture fraction 
distribution. This approach has been specifically developed for the simulation of 
turbulent diffusion flames and similar reaction processes and offers many benefits over 
the finite rate formulation approach. In the conserved scalar approach, turbulence effects 
are accounted for with the help of a probability density function or PDF. Reaction 
mechanisms are not explicitly defined by you, instead the reacting system is treated via 
either chemical equilibrium calculations or a flame sheet (mixed-is-bumed) approach. 

3. Premixed Front Tracking Formulation 

This model has been specifically developed for combustion systems, or other 
reacting systems, that are of the purely premixed type. In these problems perfectly mixed 
reactants and burned products are separated by a "flame front". The so called "G 
equation" is solved to predict the position of this front. The influence of turbulence can 
be accounted for by means of a relationship between laminar and turbulent flame speed. 

5.1 Modeling Species Transport and Finite Rate Chemistry 

FLUENT can model the mixing and transport of chemical species by solving 
conservation equations describing convection, diffusion, and reaction sources for each 
component species. This modeling capability and the inputs you provide when using it 
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are described in this section. Note that you may also want to consider modeling your 
reacting system using the mixture fraction approach. 

5.1.1 Species Transport Equations 

When you choose to solve conservation equations for chemical species. FLUENT 
predicts the local mass fraction of each species, m,„ through the solution of a convection- 
diffusion equation for the /’th species. This conservation equation takes the following 
general form: 


!■(**"' i+R ‘ +s - (A-i4) 

where R, , is the mass rate of creation or depletion by chemical reaction and S,-, is the rate 
of creation by addition from the dispersed phase. An equation of this form will be solved 
for N-l species where N is the total number of fluid phase chemical species present in the 
system. 


Mass Diffusion in Laminar Flows: , is the diffusion flux of species which 

arises due to concentration gradients. In the simplest case, this diffusion flux can be 
written as: 


J =- pD 

i .1 ’ t 




dm 


(A. 15) 


Here D, m is the diffusion coefficient for species i' in the mixture. Note that for laminar 
flows a full multi-component diffusion equation may be selected to replace Equation 
(A. 15) and that binary diffusion coefficients may be used to compute the diffusion 
coefficients in non-dilute mixtures. 


Mass Diffusion Due to Thermal Gradients: The diffusion flux, J , ■„ may 

1 

optionally be augmented by a thermal diffusion term, - D T (also called Soret 


T dx. 


diffusion). 


Mass Diffusion in Turbulent Flows: In turbulent flows, FLUENT computes the 
mass diffusion in the form: 


J 


( 



\ 


+ 


Sc, 



where S ct is the effective Schmidt number. 


M, 

PD, 


(with a default setting of 0.7). 


(A. 16) 


Treatment of Species Transport in the Energy Equation: For many 
multicomponent mixing flows, the transport of enthalpy due to species diffusion 
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( V -E*=iK K J) can have a significant effect on the enthalpy field and should not be 

neglected. In particular, when the Lewis number Le - is far from unity, when the 

k/c v 

species involved have significantly differing heat capacities, this term cannot be 
neglected. You should also consider using second order discretization schemes to reduce 
numerical errors that can be significant in the computation of this term. 

5.1.2 The Generalized Finite Rate Formulation for Reaction Modeling 

The reaction rates that appear as source terms in Equation (A. 14) are computed by 
FLUENT from Arrhenius rate expressions or by using the eddy dissipation concept due 
to Magnussen and Hjertager. Models of this type are suitable for a wide range of 
applications including laminar or turbulent reaction systems, and combustion systems 
including premixed or diffusion flames. When required, custom reaction rates can also be 
supplied through the use of user-defined subroutines. 

Reaction Rate Calculations: The source of chemical species i' due to reaction. 
R, „ is computed as the sum of the reaction sources over the k reactions that the species 
may participate in: 

A =5X* (A. 17) 

k 

where R; ■,* is the rate of creation/destruction of species V in reaction k. Reaction may 
occur in the continuous phase between continuous phase species only, or at surfaces 
resulting in the surface deposition or evolution of a chemical species. The reaction rate, 
R, t, is controlled either by an Arrhenius kinetic rate expression or by the mixing of the 
turbulent eddies containing fluctuating species concentrations. You can also enter 
alternate reaction rates using user subroutines. 

The Arrhenius Rate: The Arrhenius reaction rate is computed as: 

r u A k n c Y ex p (- /rt ) (a. is) 

j rear tan ts 

Where v, *, = molar stoichiometric coefficient for species i in reaction k 

(positive values for reactants, negative values for products) 

Mi = molecular weight of species i ' (kg/kmol) 

A = temperature exponent (dimensionless) 

A* = pre-exponential factor (consistent units) 

Cj = molar concentration of each reactant species j’ (kmol/m 3 ) 

v j’.k = exponent on the concentration of reactant j in reaction k 

Ek = activation energy for the reaction (J/kmol) 

You input values for v rk , (5 k , a*, v y k , and E k during the problem definition in FLUENT. 
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The Eddy Breakup Model: The influence of turbulence on the reaction rate is 
taken into account by employing the Magnussen and Hjertager model. In this model, the 
rate of reaction R, * is given by the smallest (i.e. limiting value) of the two expressions 
below: 

R =-v M Ap— — — — (A. 19) 

’* '* ‘ kv Ri M n 


R . =-v M AB (A. 20) 

'* *■* ' 

Where m p represents the mass fraction of any product species, P 
rriR represents the mass fraction of a particular reactant, R 
R is the reactant species giving the smallest value of R, k 
A is an empirical constant equal to 4.0 
B is an empirical constant equal to 0.5 

The eddy breakup model relates the rate of reaction to the rate of dissipation of 
the reactant and product containing eddies, k/e represents the time scale of the turbulent 
eddies following the eddy break up model of Spalding. The model is useful for the 
prediction of premixed and diffusion problems as well as for partially premixed reacting 
flows. 


Turbulent Reacting Flows: In turbulent reacting flows, FLUENT calculates the 
reaction rates from the Arrhenius expression (Equation (A. 18)) and the eddy breakup 
model (Equations (A. 19) and (A. 20)). The limiting (slowest) rate is used as the reaction 
rate and the contribution to the source terms in the species conservation and enthalpy 
equations are calculated from this reaction rate. Energy released by or required for the 
chemical reaction is accounted for in the source term of the enthalpy equation as shown 
in Equation (A. 9). 

The introduction of the kinetic term into the rate expression for turbulent flows is 
useful as it can act as a cut-off to the mixing controlled rate when chemistry is very slow. 
However, in many practical situations, the eddy breakup model describes the limiting 
rate. This is fortunate as it allows calculations to be made without knowledge of accurate 
Arrhenius rate data, the latter being hard to obtain for many industrial fuels. 

Multi-Step Reactions: For systems involving multi-step reactions, you may 
need to take special measures if you use the eddy breakup model to compute the reaction 
rate. Consider, for example, the following sequential reactions: 

CH t +—0 2 -*CO + 2H -,0 (A. 21) 

2 

co+io, ->co 2 (A. 22) 

2 2 2 
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If the product concentration determines the reaction rate (Equation (A. 20)), as is 
the case for a premixed or partially premixed flame, the reaction rate will be unphysically 
low because the concentration of C0 2 (the final product) is not included in the m p 
summation tor Reaction (A.21). The rate of breakdown of fuel is thus affected somewhat 
by the concentration of CO, which is an intermediate product with a small concentration 
in most practical cases. The resulting reaction rate is, therefore, too low. 

To remedy this problem, you should replace Reaction (A.21) with 


CH i +|o 2 ->C0 + O.OO1C0 2 +2H 1 0 (A. 23) 

When you include C0 2 as a product in Reaction (A. 23), its concentration will be 
included in the computation of the reaction rate and the resulting rate will be more 
realistic. Since the stoichiometric coefficient for C0 2 in Reaction (A. 22) is very small, 
this method should not have any adverse effect on the overall reacting system. 

It should be pointed out that the eddy dissipation model and its constants were 
originally devised for single-step representations of chemistry. The extension to multiple 
reaction steps should be made with caution. 

Stability and Convergence in Reacting Flows: Obtaining a converged solution 
in a reacting flow can be difficult for a number of reasons. First, the impact of the 
chemical reaction on the basic flow pattern may be strong, leading to a model in which 
there is strong coupling between the mass/momentum balances and the species transport 
equations. This is especially true in combustion, where the reactions lead to a large heat 
release and subsequent density changes and large accelerations in the flow. All reacting 
systems have some degree of coupling, however, when the flow properties depend on the 
species concentrations. These coupling issues are best addressed by the use of a two-step 
solution process, as described below, and by the use of underrelaxation and multiple 
sweeps of the solver. 

A second source of convergence trouble in reacting flows is the interdependence 
of the species equations on the concentrations of the other chemical species in the model. 
This issue arises because the species conservation equation for one species may be 
dominated by the reaction source term, which is a function of other chemical species 
concentrations. Because FLUENT does not simultaneously update all species 
concentrations, this interdependence can lead to convergence difficulties. This issue 
implies an increase in the number of iterations required to yield a converged solution to 
reacting flows. It also implies that simpler reaction systems converge more easily. That 
is, a one or two reaction system involving 6 or 7 chemical species converges more easily 
than a system involving 10 species in a complex reaction sequence. 
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A third convergence issue in reacting flows involves the magnitude of the reaction 
source term. When your FLUENT model involves very rapid reaction rates (i.e. much 
more rapid than the rates of convection and diffusion), the matrix solution of the species 
transport equations becomes numerically difficult. This is because the matrices to be 
solved are no longer diagonally dominant, but depend instead on large source terms that 
may drive rapid and unstable variations in the species concentrations. Such systems are 
termed “stiff" systems and are created when you define models that involve very rapid 
kinetic rates, especially when these rates describe reversible or competing reactions. Very 
rapid chemistry may be treated more successfully using the chemical equilibrium 
calculations available in the mixture fraction/PDF modeling approach. 

Two-Step Solution Procedures: Solving a reacting flow as a two-step process 
can be a practical method for reaching a stable converged solution to your FLUENT 
model. In this process, you begin by solving the flow equations only, temporarily turning 
off the species equations. When the basic flow pattern has thus been established using 
constant properties (constant concentrations), you can enable the species equations. 

Cold Flow Simulations: In combustion systems, the two-step calculation 
procedure described above is conducted by initially solving the "cold-flow", or 
unreacting flow. This cold-flow solution provides a good starting solution for the 
calculation of the combusting system. This two-step approach to combustion modeling 
can be accomplished using the following procedures: 

• Set up the problem including all species of interest in the reacting flow, but without 
the reactions defined. 

• Solve an initial solution after first turning off the product species (and, perhaps, 
enthalpy) equations. 

• Add the reaction modeling inputs, including enthalpies of formation. 

• Turn on all equations, patch an ignition source (as described below), and resolve the 
reacting flow. 

Ignition in Combustion Simulations: If you introduce fuel to an oxidant, 
spontaneous ignition does not occur unless the temperature of the mixture exceeds the 
activation energy threshold required to maintain combustion. This physical issue 
manifests itself in a FLUENT model as well and you have to supply an ignition source to 
initiate combustion. This ignition source may be a heated surface that heats the gas 
mixture above the required threshold level. Often, however, it is the equivalent of a 
spark: an initial solution state that causes combustion to proceed. This "spark" or "initial 
solution state" is supplied by you via the PATCH command: the spark is defined by 
patching a hot temperature into a region of the FLUENT model that contains a sufficient 
fuel/air mixture for ignition to occur. Often you may need to patch both the temperature 
and the fuel/oxidant/product concentrations to produce ignition in your model. The initial 
patch has no impact on the final steady-state solution-no more than the location of a 
match determines the final flow pattern of the torch that it lights. 
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6 . 


BOUNDARY CONDITIONS 


6.1 Overview of Defining Boundary Conditions in FLUENT 

Boundary Condition Options: Boundary conditions provide FLUENT with 
information on flow/thermal conditions at the boundaries of your physical model. 
FLUENT provides a wide variety of boundary condition options, including: 

• Flow inlets and exits 

• Wall boundaries 

• Symmetry boundaries 

• Periodic boundaries 

• Cyclic boundaries 

Flow inlets and exits can be defined via pressure and/or velocity specification. 
Flow exits can alternately be defined in terms of zero normal gradient (or extrapolated) 
conditions. 

Wall boundaries can be stationary or moving (e.g., rotating or sliding), slip or 
non-slip, smooth or rough. Walls may be treated via a variety of thermal boundary 
conditions (fixed temperature, fixed heat flux, or fixed external heat transfer conditions). 
Note that when an energy balance (conduction equation) is computed within wall regions, 
the thermal treatment of the wall is as a special "conducting wall" region and not as a 
boundary condition. Despite this, inputs that govern the conduction equation in walls are 
defined via the boundary condition input procedure and are covered in this chapter. 

Symmetry and cyclic boundary conditions provide a means by which the scope of 
the computational model may be reduced by exploiting the repeating nature of the 
geometry and flow pattern. Periodic boundaries allow you to model fully developed 
flows. 


Input of Boundary Conditions: Boundary conditions are defined in FLUENT by 

• assigning boundary conditions to a particular type of control volume defined 
as the "cell type" 

• assigning boundary condition values to individual control volumes via 
"patching" 

• assigning "fixed" values to individual control volumes 

In the first approach, boundary conditions are defined for all cells of a selected 
cell type. Using this approach, the boundary condition values may be uniform for all 
cells of the selected cell type, or values may vary via a polynomial function, a piecewise- 
linear spatial function, or a harmonic function. In the second approach patching- 
numerical values of boundary conditions are assigned on a cell-by-cell basis, allowing 
complete flexibility for the input of non-uniform boundary conditions. Patching can also 
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be used to input velocity boundary conditions via local coordinate systems. Fixing allows 
you to fix velocity values, for example, at live cells, simulating the effect of a fan or 
impeller without requiring the input of additional (and unknown) variable information at 
these cells. 

6.2 Flow Inlets and Exits 

FLUENT provides three types of boundary cell types for the specification of flow 
inlets and exits: velocity inlet, pressure inlet, and outlet. This section provides an 
overview and introduction to flow boundaries. 

6.2.1 Inlet and Exit Boundary Conditions for Compressible Flows 

Special attention should be paid to the setting of boundary conditions at flow 
inlets and exits when solving compressible flows. FLUENT provides several possible 
combinations of compressible flow boundary conditions. 

Supersonic Flows with supersonic Inflow and Outflow: When the flow is 
known to be supersonic at all inlets and Supersonic at the flow exit, you must choose the 
following inlet/exit boundary conditions: 

• Velocity inlet boundaries to define the velocity and static pressure (v, p s , T„ ) at 
all supersonic inlets (SUPERSONIC INFLOW BOUNDARY) 

• Outlet boundary with the command SUPERSONIC OUTFLOW 
BOUNDARY enabled to define only a single flow exit. 

Supersonic Inflow with Subsonic Outflow: When the inlet flow is supersonic 
but the flow exit is subsonic you can use the following combination: 

• Velocity inlet boundaries to define the velocity and static pressure (v, P s , To) 
at the flow inlet 

• Pressure inlet boundaries to define the exit static pressure (P s ). 

Ill-Posed Supersonic Inflow Conditions: Supersonic flow inlets should not be 
defined via input of stagnation conditions ( P n , T„) at pressure inlet boundaries. Figure 
A. 17 illustrates this ill-posed inlet treatment for supersonic flow. Such problem 
definitions do not uniquely define the static pressure (and mass flow rate) in the fluid. 


Multiple Supersonic/Subsonic Flow Inlets: Multiple supersonic flow inlets can 
be defined via (v, P 5 „ To) boundary conditions. You can also combine a supersonic flow 
inlet that has been defined via (v, P s „ To) with a subsonic flow inlet that has been defined 
via ( P„, T a ) stagnation boundary conditions or via (pv, T n ) mass flux boundary conditions. 
Note that only a single supersonic flow exit should be defined. 
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6.3 Velocity Inlet Boundary Conditions 


Velocity inlet boundary conditions are used to define the flow velocity, along 
with all relevant scalar properties of the flow, at flow inlets. In special instances, inlet 
cells may be used to define the flow velocity (but not the scalar properties) at flow exits. 


6.3.1 Inputs at Velocity Inlet Boundaries 

At a flow inlet, velocity magnitudes for all velocity components are defined. In 
addition, appropriate scalar quantities are requested as input, depending upon the problem 
scope defined. Thus your inputs may include: 

• Cartesian (or, optionally, normal, tangential, cylindrical-polar, or angular) 
velocity components 

• Mass flux (for compressible flows only) 

• Turbulence intensity/length scale 

• Temperature (total temperature for compressible flow) 

• Chemical species mass fractions (or mole fractions) 

• Static pressure (supersonic inflow boundaries only) 

7. USING THE SOLVER 


7.1 Basics of the Overall Solution Algorithm 


FLUENT solves the governing partial differential equations for the conservation 
of mass, momentum, energy and chemical species in a general form which can be written 
in Cartesian tensor notation as: 




dtp 

dx, 


+ S„ 


(A. 24) 


where <j) is the conserved quantity and the terms are the convection (LHS), diffusion, and 
source terms. The equations are reduced to their finite -difference analogs by integration 
over the computational cells into which the domain is divided. Integration in time is fully 
implicit. 

After integration of equations of the form of Equation (A. 24), the resulting 
algebraic equations can be written in the following common form: 

«VE(A-Sd=X(-V>)+S» (A. 25) 

l l 

where the summation is over the neighboring finite difference cells i = N, S, E, W, F, B 
(which stand for North, South, East, West, Front, and Back). The A’s are coefficients 
which contain contributions from the convective and diffusive fluxes and Sc and S p are 
the components of the linearized source term, S^, = S c + S p tpp. A power-law differencing 
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scheme (or the optional Second Order Upwind or Quadratic Upwind Scheme) is used for 
interpolation between grid points and to calculate the derivatives of the flow variables. 
The set ot simultaneous algebraic equations is solved by a semi-implicit iterative scheme 
which starts from arbitrary initial conditions (except at the boundaries) and converges to 
the correct solution (i.e., that which satisfies the governing equations) after performing a 
number of iterations. 

Each iteration consists of the steps, which are illustrated in Fig. A. 18. These steps 
are outlined below: 

l.The u, v, and w momentum equations are each solved in turn using current 
values for pressure, in order to update the velocity field. 

2.Since the velocities obtained in Step I may not satisfy the mass continuity 
equation locally, a "Poisson-type" equation is derived from the continuity 
equation and the linearized momentum equations. This "pressure correction" 
equation is then solved to obtain the necessary corrections to the pressure and 
velocity fields such that continuity is achieved. 

3. The K and 8 equations are solved using the updated velocity field (for turbulent 
flow only). 

4. Any auxiliary equations (e.g., enthalpy, species conservation, or any additional 
turbulence quantities) are solved using the previously updated values of the 
other variables. 

5. The fluid properties are updated. 

6. A check for convergence of the equation set is made. 

These steps are continued until the error has decreased to a required value. When 
interphase coupling is included, the source terms of the appropriate gas flow equations 
are augmented. 

7.2 Residual Reporting 

Introduction: The process of obtaining a converged solution is of great 
importance in FLUENT simulations. So that you can monitor this process, FLUENT 
provides a running report of the residuals for each equation at each iteration. The 
residuals are a measure of how closely each finite difference equation is balanced, given 
the current state of the solution. In this section, a definition of the residuals is given, and 
the use of commands for displaying un-normalized and/or peak residuals and for plotting 
residuals during the calculation is described. 
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Definition of Residuals: At each iteration of its solution algorithm, FLUENT 
reports a residual for each equation that has been solved. These residuals provide a 
measure of the degree to which each equation is satisfied throughout the flow field. 
FLUENT computes residuals for each conservation equation by summing the imbalance 
in the equation for all cells in the domain. A detailed description of the calculation of the 
residuals is provided below. 

After discretization, the conservation equation in two dimensions for a general 
variable (J) p , can be wntten as: 

Ap0 p — A E 0 E + A w (p w +A n 0 n + A s (p s + S c (A. 26) 

where 

A P =A E +A W + A y +A s -S P (A. 27) 

The nomenclature for these equations is displayed in Fig. A. 19, which shows a 
typical computational cell surrounding the node P, with neighboring nodes E, W, S and 
N. The quantity (J) might be any dependent variable (u velocity, turbulence energy, 
enthalpy, etc.). The coefficients A E , A w , A N , and A s are the finite difference coefficients 
which combine convection and diffusion through the control volume surrounding point P. 
The quantities S c and S p are components of the linearized source term which incorporate 
any terms in the equation which do not fall into the convection/diffusion form. 

The residual, R, computed by FLUENT, is the imbalance in Equation (A. 26), 
summed over all of the computational points P: 

R— ^ I A M0N +A S </) S + S c — Ap<j)p ]| (A. 28) 

nodes P 

R is the un-normalized residual and always has the SI units kg-(units of <»/sec. 

When the residuals are normalized, normalization is by the left hand side of 
Equation (A. 26), again summed over the computational nodes: 


^ nodes P 


Sib E^E “I" A n <!> n + A s </> s + S c Ap(f)p ]| 




(A. 29) 


nodes P 


In the normalized residual expressions for the momentum equations, the 


denominator term A P 0 P is replaced by A p ^ p + v p + 


w L 


R is the normalized residual. FLUENT reports normalized residuals unless you 
request that un-normalized residuals be displayed. 
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Definition of Pressure Residual: The preceding equations. Equations (A. 28) and 
(A. 29), are the residuals reported by FLUENT for all variables except the pressure. The 
residual tor pressure is actually the imbalance in the continuity (pressure correction) 
equation: 


R = X |( c w -c E + c s -c N ] 

nodes P 


(A. 30) 


where Cw, CE, Cs, CN are the convection of mass (kg/s) through each face of the control 
volume surrounding point P. Normalization of the pressure residual is accomplished by 
dividing by the residual (or continuity imbalance) at the second iteration: 


R = 


R, 


R 


iteration 2 


(A.31) 


Equation (A.31) implies that if the imbalance in continuity at iteration 2 is quite 
small, the normalized residual will be relatively large. If you make a good initial guess of 
the flow field, therefore, you can actually cause the normalized pressure residual to be 
larger than the normalized residual resulting from a poor initial guess. In such cases, you 
should examine the reduction of the un-normalized residual with iteration, since this 
value is insensitive to the initial guess. 


7.3 Judging Convergence 


Introduction: A FLUENT calculation is converged when all governing equations 
are balanced at each point in the solution domain. This section provides guidance on how 
to judge the convergence of your solution via residual values and how to monitor the 
progress towards convergence via residual histories and histories of solution variables. 

General Guidelines for Residual Values: The residuals for each flow variable 
give you a measure of the magnitude of the error in the solution at each iteration. As 
discussed in the preceding section, these residuals are normalized unless you request that 
they be un-normalized. Generally, a solution is well converged when the normalized 
residuals are on the order of lx 10° An important exception is the enthalpy residual which 
should be about 1x10 6 . (The P-1 radiation residual should also be 1 x 10 °.) 

In addition you may find that the residuals of the species transport equations will 
need to decrease to 1 x 10 5 to I x 10 6 when you are solving problems involving mixing 
of two species of very different molecular weights (e.g., H 2 and WF 6 ). If the residuals 
have decreased to this level, are monotonically decreasing, and the flow field looks 
unchanged from the solution 50 iterations earlier, then the solution can be called 
"converged". Sometimes you may not need to generate a completely converged solution, 
if you can pickup the basic features of the flow field right away. When you are interested 
in quantitative results, however, complete convergence of the solution is essential. 
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7.4 Solution of the Continuous Phase Equations 

FLUENT uses a control volume based technique to solve the conservation 
equations for mass, momentum, energy, species, and turbulence quantities described in 
the preceding sections. This control volume based technique consists of: 

• Division of the domain into discrete control volumes using a general 
curvilinear grid, 

• Integration of the governing equations on the individual control volumes to 
construct the algebraic equations for discrete unknowns (velocities, pressure, 
scalars) 

• Solution of the discretized equations 

The discretization of the differential equations, and the techniques used by 
FLUENT to solve them are described in this section. 

7.4.1 The Control Volume Technique 

FLUENT uses a control volume based technique to convert the differential 
conservation equations to algebraic equations which can be solved numerically. This 
control volume technique consists of integrating the differential equations about each 
control volume, yielding a finite-difference equation that conserves each quantity on a 
control-volume basis. 


The Non-Staggered Control Volume Storage Scheme: FLUENT defines the 
discrete control volumes using a non-staggeredgrid storage scheme as illustrated in Fig. 
A.20. In this scheme, the same control volume is employed for integration of all the 
conservation equations and all variables (pressure, Cartesian velocity components, 
Reynolds stress components, and all scalars) are stored at the control volume cell center. 


Volume Integration of the Differential Equations: The integration of the 
differential equations can be illustrated most easily in simple Cartesian coordinates, and 
is demonstrated below for a one-dimensional equation set. Consider first the one- 
dimensional differential equations for continuity, momentum, and a scalar quantity <j): 

T(/»)=0 (A. 32) 


*.(,»)=- * + i. 

dx dx Sx 
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r § 0 _ 

Sx 


+ S„ 


(A. 33) 


(A. 34) 


Equations (A. 32) through (A. 34) can be volume integrated about a control volume 
employing the Divergence Theorem: 
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(A. 35) 


f — [pu )dV = f pu . dA 

J, m ,v dx { 

Volume integration of (A. 32) on the control volume of Fig. A. 21 thus yields: 

{puA) e ~{puA) w =0 
or: 

•W„=0 

Integration of the momentum Equation (A. 33) yields: 

(pu z A) e - (pu 2 a) w =-(p t - p„ )A + 


(A. 36) 
(A. 37) 
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or: 


J e U e-J» U w =-(p e ~P w )A + 


-^-(«£ -Uph^-iUp ~“w ) 

Ax, Ax, , 


A + FAV 


(A. 39) 


and integration of the scalar equation, (A. 34), yields: 


JA-JJ. 


O _<f) 
J_E J_P 

Ax„ 


p 

Ax, , 


A+S^AV 


(A. 40) 


Note that the equations solved by FLUENT are extensions into three-dimensional 
curvilinear coordinates of those shown above for one-dimensional Cartesian coordinates. 


Closure of that Discrete Equations: the Equations (A. 37), (A. 39), and (A.40) 
are algebraic equations can be solved provided that the unknowns (u, p, and 0) are 
interpolated in a manner that relates their values at the control volume faces to the stored 
values at the control volume centers. The discretization procedures used by FLUENT to 
perform this interpolation are described in the following section. 

7.4.2 Discretization Procedures in FLUENT 

Solution of the finite-difference equations presented in the preceding section 
requires: 


• Calculation of the pressure stored at the control volume faces (p e , p w ) 

• Calculation to determine the face fluxes (J e , J w ); 
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Interpolation to relate the face values of the unknowns (u and 0) to the stored 
values at the control volume centers. 




In FLUENT, the face fluxes are obtained such that the face velocities obey an 
averaged momentum balance and the face pressures are obtained such that the velocities 
stored at cell centers obey the mass balance. These calculation procedures ensure that the 
FLUENT formulation avoids oscillatory (or "checkerboard") pressure or velocity fields 
and that the physical variation of pressure and momentum between cell centers is 
accurately represented. 

The interpolation to determine face values of the unknowns is accomplished via 
either the Power Law. blended Second Order Upwind/Central Difference, or QUICK 
interpolation schemes, described below. 


The Power Law Scheme: The power law interpolation scheme interpolates the 
face value of a variable, 0, using the exact solution to a one-dimensional convection- 
diffusion equation. This one-dimensional equation describes the flux of 0 as: 


ox ox ox 


(A. 41) 


where T and pu are constant across the interval 5x. Equation (A.41) can be integrated to 
yield the following solution describing how 0 varies with x: 


<?(*)- 

0L-0O 


exp 

f 

P t ~ -1 

I 

kl 


exp(Pj-l 


where: 


<t>„ - <j)\ t=0 and P e is the peclet number, 

<Pl=<P \„L 

p 

e r 


(A. 42) 


(A.43) 


FLUENT uses 7.4-11 in an equivalent "power law" format, as its default 
interpolation scheme. 

The variation of 0(x) between x = 0 and x = L is depicted in Fig. A. 22 for a range 
of values of the Peclet number. Figure A. 22 shows that for large Pe, the value of 0 at x = 
L/2 is approximately equal to the upstream value. This implies that when the flow is 
dominated by convection, interpolation can be accomplished by simply letting the face 
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value of a variable be set equal to its "upwind" or upstream value. When the Pe = 0 (no 
flow, or pure diffusion), Figure A. 22 shows that 0 may be interpolated via a simple linear 
average between the values at x = 0 and x = L. When the Peclet number has an 
intermediate value, the interpolated value for 0 at x = L/2 must be derived by applying 
the "power law” equivalent of Equation (A. 42). 

Higher Order Interpolation Schemes: As an alternative to the Power-Law 
Differencing Scheme, FLUENT provides two higher order schemes: QUICK and a 
blended Second Order Upwind/Central Difference scheme. These schemes compute the 
face value of an unknown (e.g., 0 f ) based on the values stored at the two adjacent cell 
centers (0p and 0 E ) and on a third cell center at an additional upstream point (e.g., 0 W ). 
Using the nomenclature depicted in Fig. A. 23, the face value can be written in terms of 
these neighbor values as: 
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(A.44) 


where 0f is the face value, 0 D is the downstream value, 0c is the center cell value, and 0c 
is the upstream value. Equation (A. 44) computes the face value with second or third 
order accuracy, depending upon the choice of 0, as shown in Table (A.l). FLUENT 
chooses 0 for the two higher order schemes in a manner which eliminates oscillations or 
overshoots, as noted below. 

Bounding of the Higher Order Schemes: While higher order schemes provide 
greater accuracy, numerical instabilities can occur unless the interpolation is 
appropriately bounded. Second-Order Central Difference will yield oscillations at 
discontinuities, and Second-Order Upwind or QUICK will produce undershoots and 
overshoots. Linear combinations of the methods can, however, produce a scheme which 
is both stable and accurate. 

Bounding of the Second-Order Scheme: For the Second-Order scheme, such 
linear combinations are created through variation of the parameter 0 in Equation (A.44). 
In the blended Second-Order Upwind/Central Difference Scheme used by FLUENT, 0 is 
computed as: 

6> = max(o,min(0 c , l)) (A.45) 

where 


0c = 


0c 0u 

0D~0U 


(A.46) 
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Bounding of the QUICK Scheme with Limiters: For a uniform grid. Equation 
7.4-13 can be rewritten as 

<P, =<Pc+\[{i + k){<pD-<t>c) + {i-k){<Pc-0v)] (A. 47) 

where the numerical parameter k controls the order of the scheme. Unbounded terms 
corresponding to central differencing (CD), second order (linear) upwind differencing 
(LUDS), and QUICK schemes arise by setting k = 1, -1, and 0.5, respectively. To 
provide a monotonic scheme in terms of Equation (A.47), a slope limiter (r f ) is 
introduced: 


=«V +^[(l + ^V( r /)(^D -<t>c)+ O + - <Pu )] 


(A. 48) 


where 
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r, = — %s- 
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(A.49) 


0 C is as defined in Equation 7.4-15. 

In FLUENT, for the QUICK scheme (k = 0.5), four different kinds of limiters 
have been implemented: 

• UMIST or SMART: 


^(r)=max[0,min(2r , 0.75 + 0.25r , 2)] 


• MUSCL: 


V{r)= i 
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compressible MINMOD: 


yr(r)=max[0,min(l ,5r)] 


• SUPERBEE: 


y/{r) = max [0 , min(2r , l), min(r,2)] 


(A. 50) 


(A. 51) 


(A. 52) 


(A. 53) 
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These limiters provide monotonic behavior of the variables in regions of steep 
gradients, and can be characterized by their different convergence properties: the steeper 
the limiter, the more difficult it is to converge. UM1ST or SMART (the default) has 
optimal characteristics in terms of convergence and boundedness, while SUPERBEE 
(which is very good at resolving discontinuities) is the most difficult to converge. 
FLUENT also allows you to use the unbounded form of QUICK w'ithout limiters (v|/(r) = 
1 ). 


Linear Pressure Interpolation: For supersonic compressible flows in which 
shocks can occur, FLUENT provides a linear pressure interpolation option. The standard 
pressure interpolation at cell faces is a physical interpolation, which works well for nearly 
all flows, but which may fail to satisfy conservation of total temperature in compressible 
flows with shocks. Note that linear pressure interpolation cannot be used with porous- 
media flows or flows with strong body forces. 

Higher-Order Density Interpolation: Four interpolation schemes are available 
for the calculation of density: upwind (default), linear, second-order-upwind, and 
QUICK. The default upwind scheme sets the density at the cell face to be the upstream 
cell-center value. This scheme provides stability for the discretization of the pressure- 
correction equation, and gives good results for most classes of flows. The linear density 
interpolation scheme uses a simple linear interpolation. For some incompressible flows 
(e.g., natural-convection flows or flows with large changes in density due to temperature 
or species concentration), better results can be obtained by using the linear density 
interpolation. In compressible flow with shocks, the default upwind scheme may tend to 
smooth the shocks; you should use the second-order-upwind or QUICK scheme for such 
flows. For compressible flows with shocks, using the QUICK scheme for all variables, 
including density, is highly recommended. 

Viscosity-Weighted Velocity Interpolation: To calculate the differential viscous 
terms for the momentum equations and for the generation of turbulence, FLUENT 
provides an optional viscosity-weighted velocity interpolation scheme in addition to the 
default linear interpolation scheme. The viscosity-weighted method is based on the 
continuity of shear stress for control volumes. For the simple ID case shown in Fig. 
A. 24, the velocity at face c obtained with the viscosity-weighted interpolation will be 


(Ax) t _ 

Mp{ A*) e+ + Me (A* )<r- 


(A.54) 


The viscosity-weighted interpolation option will give better results for turbulent 
kinetic energy inside shear layers without excessive generation of turbulence. 
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7.4.3 The Pressure-Velocity Coupling Algorithms: SIMPLE and SIMPLEC 

The continuity and momentum equations, in three dimensions, provide four 
equations for solution of four unknowns: u,, u 2 , u 3 , and p. Simultaneous solution of this 
(linearized) equation set would provide a solution in which all four unknowns satisfy 
each of the (linearized) mass and momentum equations. Because a simultaneous solution 
is computationally intensive, in FLUENT the equations are solved sequentially. In the 
sequential solution process, an equation describing the update of pressure is required, and 
is not explicitly available via the mass or momentum balances. The SIMPLE family of 
algorithms is based on using a relationship between velocity and pressure corrections in 
order to recast the continuity equation in terms of a pressure correction calculation. 

Coefficient Form of Momentum Equation: the The discretized one-dimensional 
momentum equation developed in Section 7.4.1 above, can be written in a shortened 
notation as 

A p u p =^A, s h vs + ( Pw - p e )A + S (A. 55) 

MB 

where Ap and A NB are coefficients containing the convection and diffusion contributions 
in the momentum equation. Equation (A. 39), and the subscript NB refers to neighbor 
points (e.g., E, IT in ID or E, W, N, S in 2D). 

The SIMPLE Algorithm: The SIMPLE algorithm starts with substitution of a 
guessed pres sure field, p*, into the momentum equations (Equation (A. 55)) which can 
then be solved to obtain a "guessed" velocity field, u*: 

A P u p =X A m» m 1 + ip\ -P*) A + 5 (A. 56) 

MB 


Equation (A. 55) is solved for the "guessed" velocity field (e.g. w* ). The actual velocity 
and pressure fields are related to the "guessed" values, u p and p*, and 

u p =Up+u p (A. 57) 

Pe=P'e+Pl (A. 58) 

where w* and p* e are the velocity correction and pressure correction, respectively. 

Substitution of Equations (A.57) and (A.58) into Equation (A.55), followed by 
subtraction of (A. 56), yields a “momentum balance" in terms of the velocity and pressure 
corrections: 

A p Up — ^ A nb u nr + [p w ~ p , )a (A. 59) 

MB 

Equation (A. 59) is used to relate the pressure and velocity corrections as 
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(A. 60) 


u 'p =-^-(pi - p'M 

A p 

where the term containing the neighbor influences (Xvs A nb u 'nb) is simply dropped for 

convenience (and will be zero at convergence since it involves only velocity corrections 
at the neighbor points) . 

Equations similar to Equation (A. 60) are used to cast the continuity equation in 
terms of an equation for the pressure correction. The basic mass balance 

(piiA) e -{piiA) K =Q (A. 61) 

is first written in terms of the velocity u* + u’: 

{pA) e (u’+u'l- [pA) w (u * + u ' )„. = 0 (6.62) 

Using equations like Equation (A. 60), the continuity equation can be recast as an 
equation for the pressure correction: 

(pAu r l - (pAu * ) w + {pA\ (p'r-p'e)- {pA) w y - (p' w ~ p' P )=0 ( A.63) 

Equation 7.4-32 can now be solved for a correction to the pressure field which is 
then used to compute the velocity correction via Equation (A.60). Finally, the velocity 
and pressure are updated via Equations (A. 57) and (A. 58). 

The SIMPLEC Algorithm: The SIMPLEC algorithm is a variant on the standard 
SIMPLE algorithm described above. The derivation of SIMPLEC is the same as that for 
SIMPLE from Equations (A. 55) - (A. 59). From that point on, the SIMPLEC derivation 
proceeds as follows. 

Equation (A. 59) is used to relate the velocity and pressure corrections, after 
subtraction of a new term, A NB u' NB ), from both the left and right sides: 



(A. 64) 


Next, the term involving the difference [u' NB -u p ) is dropped, under the argument 

that this term is small and will vanish at convergence when the corrections are zero. The 
resulting relationship between pressure and velocity correction becomes 
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(A. 65) 


u D =- 
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Equation (A. 65) can be contrasted to the correction relationship used in SIMPLE 
(Equation (A. 60)) by noting that here the neighbor term which is dropped involves a 
difference term(^ A NB [ii‘ NB -u' p )) which is small compared to the neighbor term 

(X vfl Avb ) dropped in SIMPLE. 


Equations similar to Equation (A. 65) are used to cast the continuity equation in 
terms of an equation for the pressure correction. The mass balance equation (in ID 
Cartesian form) is first written in terms of the velocity h* +u' : 

(pA) e (u'+u'\ ~{pA) w («* +u' ) H , =0 (A. 66 ) 

Using equations similar to Equation (A. 65), this continuity equation can be recast in 
terms of pressure correction as 
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(A.67) 


Equation (A.67) can now be solved for a correction to the pressure field which is 
then used to compute the velocity correction via Equation (A. 65). Finally, using 
Equations (A. 57) and (A. 58), the velocity and pressure corrections thus obtained are used 
to update the current velocity and pressure fields. 


7.4.4 The Iterative Solution Procedure 

The SIMPLE and SIMPLEC algorithms described above relate the velocity and 
pressure fields which satisfy the linearized momentum and continuity equations at a 
point. Because FLUENT does not solve the equations at all points simultaneously, and 
because the equations are coupled and non-linear, an iterative solution procedure is 
required with iterations continuing until all equations are satisfied at all points. 

Each iteration of FLUENT’S solution procedure consists of the following steps: 

1. The U|, U 2 , and U 3 momentum equations are solved in turn using the guessed 
pressure field, p*. 

2. The pressure correction equation (mass balance) is then solved to obtain the 
necessary corrections to the pressure field. Corresponding adjustments to the 
velocity components are also made. 
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3. For turbulent flows, the k and e equations (or Reynolds stress transport 
equations) are solved using the updated velocity field to obtain the distribution of 
the effective viscosity and/or Reynolds stresses. 

4. Any auxiliary equations (e.g. enthalpy, species conservation, and/or radiation) 
are solved using the previously updated values of other variables. 

5. Fluid properties are updated 

6. When interphase coupling is to be included, the source terms in the appropriate 
continuous phase equations may be updated via a dispersed phase trajectory 
calculation. 

These steps can be continued until the error in each conservation equation within 
each volume and hence over the global domain has decreased to a required value. 

7.4.5 Iterative Full-Field Solution of the Equations 

On those computer system supporting the required matrix libraries, FLUENT 
optionally uses a full-field iterative solver called GMRES (Generalized Minimum 
Residual Method) to solve the system of equations. This full-field solver can vield 
significant speed-up of convergence compared to the Line Gauss-Seidel technique 
described below. 

7.4.6 Line-by-Line Solution of the Equations 

The algebraic equation to be solved by FLUENT for any variable at point P may 
be written as: 


^NB^NB ^0 (A. 68) 

,VB 

where the subscript NB denotes neighbor values, the coefficients A p and A N b contain 
convection and diffusion coefficients, and is the source of <J) in the control volume 
surrounding point P. For each unknown, <]), an equation of this form must be solved at all 
points within the domain. This solution process may be accomplished via a "line-by- 
line" solver, in which the equations along a single "line" of cells (e.g. a line of I = 
constant) are solved simultaneously. The line-by-line solver gives rise to a tri-diagonal 
matrix which is solved via Gaussian elimination to update the values of <|> along the line 
considered. This procedure, also referred to as Line Gauss-Seidel (LGS), is repeated for 
all lines in the domain so that <j> is updated at all points P. FLUENT provides flexibility in 
its line-by-line solution procedure, giving the user control over the directions) of the lines 
to be considered and over the number of times each line is visited in order to update a 
given variable within each global iteration loop. 
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Multigrid and ID Block Correction to Accelerate Line-By-Line Solver: While 
the LGS solver is good at reducing local errors, it is relatively poor at reducing errors 
which require collective change of many cell the values (long-wavelength errors). 
Therefore the speed of solution of the LGS solver deteriorates with increasing grid size. 
Also for large aspect-ratio cells in combination with gradient boundary conditions, as 
may occur in highly stretched grids, convergence may be poor. FLUENT provides two 
devices to speed up the LGS procedure: 

• One-dimensional block correction 

• Multigrid (MG) acceleration. 

Both methods accelerate the convergence of the equations by deriving global 
corrections that drive the solution toward global conservation. One-dimensional block 
correction is also useful in combination with multigrid acceleration, when cells are highly 
stretched or otherwise strongly anisotropic coefficients are to be expected. Each of these 
methods is described in the following sections. 

ID Block Correction vs. Multigrid: One-dimensional block corrections reduce 
long-wavelength errors in the direction in which they are applied, but may introduce large 
short-wavelength errors. Therefore, application of block correction is followed by 
application of LGS to reduce the short-wavelength errors. Sufficient sweeps must be 
specified to this purpose, otherwise the effect of block correction may be adverse. If very 
steep gradients exist in the final solution, ID block correction may increase the solution 
effort for this reason. These disadvantages are to a large extent avoided in multigrid 
(MG) acceleration. 


7.4.7 Multigrid Acceleration of the Line-by-Line Solver 

The line-by-line solution approach is accelerated in FLUENT using an Additive- 
Correction Multigrid procedure. This procedure computes corrections to the current 
solution field within successively coarser blocks of control volumes. The corrections 
obtained on these multiple coarse-grid levels are used to refine the fine grid solution and 
to thereby accelerate convergence. 

The Additive Correction Multigrid Scheme: Additive Correction is a technique 
which seeks to apply a correction, t> / , to the current solution field, <)>*, such that the 
resulting corrected solution, 0 + ()/ , obeys a global conservation of 0 within some 
subregion considered. In FLUENT, the subregions considered are constructed via 
grouping together of neighboring control volumes. The equations to be solved on each 
multigrid level are constructed via summation of the equations on the original fine grid 
level. This summation, over the fine grid cells used to construct the coarse grid element, 
can be written as: 
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XXEk* =Am +■ a; + a;I <*w 
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+ Ay< ^,y-U + A'* #,*„] +A^ 0y*-l + A* J (A ' 69) 


where the coefficients A £ A W ,A' V , A s , A u , and A D are the neighbor coefficients (Anb) of 
the point P at cell ijk of the fine gnd. 


Next, the unknowns, 0 ]Jk , are assumed to consist of the current solution value in 
each cell, <p ' jk , and a correction, <t>‘ IJK on the IJKlh block of the coarse gnd level: 


A/* =ti jk (A. 70) 

Substitution of Equation (A. 70) into (A. 69) yields an equation to be solved for the 
correction field, <j>' lJK , on the coarse grid level. This correction equation relates the 

correction in the IJKlh block of the coarse grid to the corrections in the neighbor blocks 
on the coarse grid. The set of equations derived via similar summations on each coarse 
grid block are then solved iteratively, using an alternating direction Line-Gauss-Seidel 
solution method. The set of corrections thus obtained are then added to the current fine 
grid solution as in Equation (A. 70). The resulting corrected solution then satisfies global 
balances on the coarse grid level. Satisfaction of local conservation on the fine grid level 
then proceeds as usual, using the line solver to update <j>* k . 

Control of Grid Levels During the Multigrid Solution Process: The multigid 
procedure invokes calculations on the next coarser grid level when the error reduction 
rate on the current level is insufficient, as defined by the parameter P: 

(A. 71) 

Here Rj is the absolute sum of residuals computed on the current grid level after 
the ith sweep of the line solver on this level. Equation (A.71) states that if the residual 
(error) present in the iterative solution after i sweeps is greater than some fraction, p, of 
the residual present after the (i-l)th sweep, the next coarser grid level should be visited. 
Modifications to P control the frequency with which higher grid levels are visited. 

Provided that the residual reduction rate is sufficiently rapid, FLUENT will 
converge the correction equations on the current grid level and apply the resulting 
corrections to the solution field on the next finer grid level. The correction equations on 
the current grid level are considered sufficiently converged when the error in the 
correction solution reduces to some fraction, a, of the original error on this grid level: 

R, < a R„ (A. 72) 
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Here, R, is the residual on the current grid level after the ith iteration on this level, 
and R„ is the residual which was initially obtained on this grid level at the current global 
iteration. Equation (A. 72) is also used to terminate calculations on the lowest (finest) 
grid level during the multigrid procedure. Thus, FLUENT will continue sweeping the line 
solver on each grid level (including the finest grid level) until the criteria of Equation 
(A. 72) is obeyed (or until a maximum number of sweeps has been completed, in the case 
that the criteria of Equation (A. 72) is never achieved.) 

Adding ID Block Correction with Multigrid: As noted above, the convergence 
of the coarse grid correction equations may be accelerated by computing corrections on 
an even coarser grid and/or by applying ID block correction to the current grid level. ID 
block correction can be understood as a simple version of the Additive Correction 
Multigrid technique, with the coarser grid level consisting of a one-dimensional grid 
created by summation over all cells in two of the three grid directions on the finer grid. 
ID block correction is of particular benefit in problems with large cell aspect ratios and 
in problems with large abrupt changes in material properties. 

7.4.8 ID Block Correction 

FLUENT optionally applies a ID block correction to the continuous phase 
transport equations in order to accelerate convergence of the line-by-line solution process 
described above. In the block correction procedure used by FLU ENT, the transport 
equations for a variable 0 are summed in each plane of constant 1, J, or K, and the 
resulting set of "one-dimensional" equations are solved for a correction to on each plane. 
The details of this process are described below. 

Consider the discretized equation for <p ijk at point P in 3D: 


Ay* 0,jk + A y* 0|+1 ]k + A Ok 0.-1 Jk + Ay* 0y-+u + Ay* 0/y-i* + Ay* 0,y*+i + Ay* 0<y*-i Ay* (A. 73) 


where the superscripts E, W, N, S, U, and D represent the east, west, north, south, up, and 
down neighbors of the point P. For the case of 3D to ID block correction in the I 
direction, the equations are summed in each I plane: 

XX (a; <#,» + K *,.» + A,i +A", +A° *<*.,)+ 

' * (A. 74) 

X X a i 0,+i jk + X X aJ 0,-iy * = X X a* 

j k j k j k 


If 0 iJk is the current best estimate of the solution, the improved solution will be 
given (p l]k ~0 IJk +S,. Here 8 ( is the correction over the plane. Substituting this into the 
above equation yields 
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where: 




Ajk fyjk $i+\jk Ajkfii i jk A } k 0tj+\k Ajkfiij-u Ajkftijk+\ Ajk$ijk-\ 

J k 

) (A. 76) 

— 

Defining a set of coefficients for these ID equations, the correction equation can 
be written as 


A, p 8, +A, e S l+l +A* S -b, 

(A. 77) 

— 

where: 



a; -£E(a; +a;/» + aJ/ s + 4 /" + a“/“ ) 

J k 

(A. 78) 


a £ -sia; 

j k 

(A. 79) 


a" - six 

J k 

(A. 80) 


The factors /,/,/, and f in Equation (A. 77) are used to remove the influence of 
neighbors in the J or K directions in the I direction block correction equation. This is 
accomplished by setting/' = 1 if A is a computational cell and 0 otherwise, with similar 
treatments for/, /, and/. 

The correction equations. Equation (A. 77), for the corrections5[, at each / plane 
are solved using the Thomas tri-diagonal matrix algorithm, and then the unknowns are 

corrected using </> ijk =<j> ijk + 8, ■ Similar equations are derived and solved for the 

corrections in the J and K directions. 
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7.5 Speeding Convergence 


Introduction: Convergence can be hindered by a number of factors. Large numbers of 
computational cells, overly conservative underrelaxation factors, and complex flow physics are 
often the main causes. In this section, some of the numerical controls and modeling techniques 
that can be exercised to enhance convergence are examined. 

7.5.1 Initial Guesses and Step-by-Step Solution Processes 

Introduction: Two of the simplest techniques to get your calculation off to a good start 
are covered in this section. These include starting with guessed values for some of the variables 
at some or all locations within the domain, and techniques for solving difficult problems in 
stages. 


Patching Initial Guesses: Supplying an initial guess for important flow variables, using 
the PATCH command, is an excellent way to begin the solution of a problem. In the absence of 
an initial guess, all of the variables are assumed to have a value of 0 throughout the domain, with 
the exception of temperature, (which has the default value of the freezing point for water, 273 K) 
and the turbulence parameters k and c (which are set to averages of the inlet values). Use of the 
PATCH option is illustrated below, where the entire domain (except protected boundaries) is 
covered with the patched values for temperature (TE) and u velocity (UV) by using the defaults 
on the minimum and maximum I and J indices. 

Some of the most common examples where patching can be helpful are listed below: 

• An initial guess for the fluid temperature should always be made in problems 
involving heat transfer and/or when the gas law is used. 

• An initial guess of the velocity field is particularly useful when a large number of 
computational cells is used in the problem setup. 

• In natural convection problems, it is helpful to patch both the temperature and 
velocities at the start of the problem. 

• In problems with reacting flows, you can choose to solve the reacting flow right from 
the start or use a step-by-step method as described in the next section. If you decide to 
use the first option, it is a good idea to patch mass fractions for the reactants in the 
region where the reaction is most likely to occur. In combustion problems, where a 
threshold temperature is required to sustain a reaction, it is a good idea to patch in this 
temperature at the start. 

Patching the Heat Capacity When Patching Temperature: When you supply an initial 
guess for the fluid temperature, FLUENT will update the enthalpy using the current stored value 
for the mixture heat capacity c p . When the heat capacity is composition dependent, FLUENT’S 
current stored value for c p , may be considerably different from that expected in the converged 
solution. In some instances, this may lead to poor initialization of enthalpy, resulting in a sudden 
change of temperature, away from your patched values, when the solution process updates the 
stored value of c p . You can work around this start-up issue by first patching c p , and then patching 
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temperature. This is not required, however, as the converged solution will be independent of your 
patched initial guess. 

Step-by-Step Techniques: One important technique for speeding convergence for 

complex problems is to tackle the problem one step at a time. When modeling a problem with 
heat transfer, you can begin with the calculation of the isothermal flow. When modeling a 
reacting flow, you can begin by computing a partially converged solution to the non-reacting 
flow, possibly including the species mixing. When modeling a dispersed phase, such as fuel 
evaporating from droplets, it is a good idea to solve the gas phase flow field first. Such solutions 
generally serve as a good starting point for the calculation of the more complex problems. This 
option is explained below. 

FLUENT automatically solves each equation that is turned on in the problem setup via 
the Models panel or the DEFINE-MODELS menu. If you specify that the flow is turbulent, 
equations for conservation of turbulence quantities are turned on. If you specify that FLUENT 
should calculate temperature, the enthalpy equation is activated. Convergence can be speeded by 
focusing the computational effort on the equations of primary importance. The SELECT- 
VARIABLES table allows you to turn individual equations on or off temporarily. 

A typical example of using the SELECT-VARIABLES command is in the computation of 
a flow with heat transfer. Initially, you would define the full problem scope, including the 
thermal boundary conditions and temperature-dependent flow properties. Following the problem 
setup, you would use the SELECT-VARIABLES command to turn off the enthalpy equation. 
You could then compute an isothermal flow field, using a patched value for the temperature of 
the fluid. When the isothermal flow was reasonably well converged, you would turn the enthalpy 
equation back on. You could, in fact, turn off the momentum and continuity equations while the 
enthalpy field was being computed. When the enthalpy field began to converge well, you would 
turn the momentum and continuity equations back on so that the flow pattern could adjust to the 
new temperature field. The temperature would couple back into the flow solution via its impact 
on fluid properties such as density and viscosity. The temperature field would have no effect on 
the flow field if the fluid properties (e.g., density, viscosity) did not vary with temperature. 

Turning Reactions On and Off: To solve a species mixing problem prior to solving a 
reacting flow, you should set up the problem including all of the reaction information, and save 
the complete case file. To turn off the reaction so that only the species mixing problem can be 
solved, you can open the Define Species panel and decrease the number of Reactions to 0. 

Once the species mixing problem has partially converged, you can return back and restore 
the correct number of reactions as originally specified. You can then resume the calculation 
starting from the partially converged data. 

As noted above, for combustion problems you may want to patch a hot temperature in the 
vicinity of the anticipated reactions before you restart the calculation. 
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Starting from a Previous Solution: In the above examples, it was shown how you can 
solve a complicated problem in steps. In doing so, you can use the partially converged data from 
the first step, where some of the equations or models are not included, as a starting point for a 
later step in which these equations or models are included. This is a general technique w hich can 
be used in other ways. Suppose, for example, that you want to compare the flow through a 
channel with and without an internal obstruction. You can solve the problem first without the 
obstruction and obtain a converged solution. You can then modify your case file so as to include 
the obstruction, and read in the data from the first calculation as a starting point for the second. 

7.5.2 Underrelaxation 

Introduction: Because of the nonlinearity of the equation set being solved by FLUENT, 
it s not generally possible to obtain a solution by fully substituting the "improved" values for 
each variable which have been generated by the approximate solution of the finite difference 
equation. Convergence can be achieved, however, by underrelaxation which reduces the change 
in each variable produced during each iteration. In a simple form, the new value of the variable 
<D P at node P depends upon the old value, <Dp,oid> the computed change in <D P , A(D P , and the 
underrelaxation factor, a, as follows: 

0 p — 0 p 0 ij + ccA(/) p (A. 81) 

In FLUENT, the default underrelaxation parameters for all variables except the 
velocities are set to low values in order to ensure convergence in the largest possible number of 
cases. Unfortunately, this may not give rise to the fastest rate of convergence, and an 
improvement can often be obtained by a judicious increase in one or more of these parameters. 

Increasing the Underrelaxation Parameters: As mentioned earlier, FLUENT currently 
incorporates two different velocity-pressure coupling algorithms: SIMPLE and SIMPLEC, both 
of which have been modified to handle the complexities of nonorthogonal body-fitted 
coordinates. The optimal choice of underrelaxation factors is dependent not only on whether 
SIMPLE or SIMPLEC is used, but also on the type of flow being modeled and properties of the 
mesh. Underrelaxation factors for the transport equations are set in the range of 0 to 1. The 
smaller the factor, the heavier is the degree of underrelaxation and the greater is the degree of 
control exercised over the change permitted from one iteration to the next. While small 
underrelaxation factors damp out nonlinearities, they also serve to inhibit the rate of convergence 
for more straightforward problems. 

The default underrelaxation factors in FLUENT are set rather conservatively. In other 
words, the low default values ensure stability for virtually all problems, but in many cases lead to 
slow convergence. Provided that the coupling between the transport equations is reasonable, 
property variations are not very large, and the mesh is not very distorted, larger values of the 
underrelaxation factors will lead to faster convergence. Table A.2 below summarizes suggested 
underrelaxation factors for problems that are reasonably benign and for time-dependent flows 
where the advancement of the solution from one time step to another already provides a 
reasonably good guess. 
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Note that the only difference between the SIMPLE and SIMPLEC settings is the 
underrelaxation factor for the pressure correction equation. When SIMPLEC is used, generally 
no underrelaxation is required for the pressure correction equation (i.e., the factor should be set at 
1.0). An exception to this may be required for highly compressible flows where the effect of 
large density variations on stability can be controlled by decreasing the underrelaxation factor for 
the pressure correction equation. 

The residuals should decrease monotonically or near monotonically after an initial startup 
phase where some residual increase or oscillations may be observed (typically 10-20 iterations). 
When the problem being solved is particularly "bard", however, one or more of the residuals will 
not decrease and corrective measures are required. These "hard" problems probably have one or 
more of the following properties: 

• mesh distortion 

• large body forces (e.g., due to buoyancy, rotation) 

• multiphase phenomena (e.g., large particle loading for the Lagrangian model, Eulenan 
multiphase model) 

• large property variations (e.g., in combustion or phase change problems with highly 
temperature- and composition-dependent properties, non-Newtonian flows) highly 
compressible flows 

The corrective measures usually revolve around reduction of the underrelaxation factors. 
Reduction should be done in increments of 0.05 to 0.1 and applied to the equations for which the 
residuals are increasing. If the offending residuals are for the pressure and momentum equations, 
the reduction should be applied to the velocities. Increasing the number of sweeps may also be 
worthwhile. For turbulent flow problems, combustion problems, and problems with variable 
properties, the underrelaxation factors for enthalpy, species, and viscosity can be reduced as well. 

Occasionally, you may make changes in the underrelaxation factors and resume your 
calculation, only to find that the residuals begin to increase. This often results from increasing the 
underrelaxation factors too much. A cautious approach is to save a Data File before making any 
changes to the underrelaxation factors, and to give the solution time to adjust to the new 
parameters. Typically, an increase in the underrelaxation factors brings about a slight increase in 
the residuals, but these increases usually disappear as the solution progresses. If the residuals 
jump by a few orders of magnitude, you should consider halting the calculation and returning to 
the last good Data File saved. 

7.5.3 Solution Methods for the Discretized Equations 

FLUENT uses iterative techniques to obtain an approximate solution to the set of 
discretized equations. Iterative methods of solution require much less effort than a direct (exact) 
solution, but attention has to be paid to solution accuracy at each iterative stage. Controlling 
solution accuracy requires some basic understanding of equation solution procedures available in 
FLUENT. 


48 



The Line-by-Line Solver (LGS): One technique used by FLUENT is the line-by-line 
solution tech nique, known as Line-Gauss-Seidei (LGS). In LGS the equations are solved 
simultaneously for small groups of cells, one at a time. The groups consist of lines of cells, w here 
a line is either a complete row or a complete column of cells. The direction of the line is called 
the sweep direction. During the update of the solution on a line the solution on neighboring lines 
is treated as correct and kept constant. The direction in which the solution process proceeds 
through the domain (line-by-line in 2D, or plane-by-plane in 3D) is called the marching direction. 
A solver sweep is a complete traversal of the domain. A number of sweeps may be required for 
sufficient solution accuracy. 

Multigrid and ID Block Correction to Accelerate the Line-By-Line Solver: Generally 
LGS is good at reducing local errors, but poor at reducing errors which require collective change 
of many cell values (long-wavelength errors). Therefore the speed of solution of LGS deteriorates 
with increasing grid size. Also for large-aspect-ratio cells in combination with gradient boundary 
conditions, as may occur in highly stretched grids, convergence may be poor. FLUENT provides 
two devices to speed up the LGS procedure: 

• One-dimensional block correction. 

• Multigrid (MG) acceleration. 

Of these, MG acceleration is recommended for pressure in almost any case and is active, 
by default, on the pressure and enthalpy equations. 

GMRES: A Full-Field Iterative Solver: On selected computer hardware, FLUENT also 
provides access to a full-field iterative solver, GMRES. This alternative to the line-by-line solver 
is provided on those computers which provide access to GMRES math libraries. 

7.5.4 Controlling the Sweeps of the Line-by-Line (LGS) Solver 

Introduction: Sometimes a problem shows extremely slow convergence, with the 
residuals for one or more equations remaining nearly constant over a large number of iterations. 
In addition to modifications in the underrelaxation parameters, two other types of modifications 
in the solver parameters are available which can help: 

• The sweep direction, which controls the way in which cells are grouped together to be 
solved 

• The number of sweeps on each equation, which controls the extent to which the 
solution of any single equation is refined during each iteration 

In this section, the nature of these solution controls is explained in more detail, along with 
guidelines on how and when to make changes to them. 

Solver Sweep Direction: If an exact solution were performed on a computational domain 
with N cells, an N x N matrix would need to be solved. Rather than solve a problem via this 
"direct method," FLUENT uses an iterative solution technique, solving the equations 
simultaneously for small groups of cells, one group at a time. The groups contain either a 
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complete row or a, complete column of cells, starting at one boundary and ending at the other. 
The solver sweep direction is the direction in which these groups of cells are formed, i.e., 
columns represent one direction while rows represent another. To further illustrate the meaning 
of the solver sweep direction, consider the computational grid shown in Fig. A. 25. 

Figure A. 25 illustrates a solver sweep in the J-direction (direction 2). Sweep Direction of 
J means that the equations are solved simultaneously for a single column of cells (or a sweep of 
cells in the J-direction). The solver then marches forward in the positive Indirection to solve the 
next column of cells. 

Choosing the Best Sweep Direction: Generally speaking, the sweep direction should be 
in a direction that is normal to the primary direction of flow. In this manner, the effects of the 
boundary conditions are propagated along the sweep direction, and this information can then be 
propagated forward along the direction of flow. Thus the choice of sweep direction 2 (or J) works 
well when flow is primarily in the I-direction, as indicated in Fig. A. 25. When flow is primarily 
in the J-direction, you should consider Sweep Direction of I, as shown in Fig. A. 26. 


Alternating Sweep Direction: FLUENT uses an alternating sweep direction, by default, 
sweeping lines in the I direction and then sweeping lines in the J direction (in 2D). 

The alternating-direction solver begins by using the first sweep direction to update the 
current dependent variable at each point. It then repeats the calculation using the alternate 
direction as the sweep direction. The process is illustrated in Fig. A. 27. This technique is 
particularly useful for problems in which there is no single dominant flow direction, e.g., flows 
with jets or particles injected into a crossflow. 

If the option is used in 3D, the solver operates as in the following example. Suppose you 
choose MARCHING DIRECTION 3 (March Direction of K). The alternating-direction solver 
then updates the dependent variables by using: 

1. sweep direction 1(1), marching direction 3 (K) (first) 

2. sweep direction 2 (J), marching direction 3 (K) (second) 

In other words, the marching direction is always the same, and the sweep direction 
alternates (beginning with the smaller directional index). 

Sweep and March Directions in 3D Problems: In 3D problems, you can choose both a 
sweep direction and a marching direction. The sweep direction is as defined above, and the 
marching direction is that along which the sweep solution propagates. The rule of thumb once 
again is to march in the flow direction and sweep in the direction across which boundary 
conditions are imposed. For example, SWEEP DIRECTION 3 (Sweep Direction of K) should be 
selected when boundary conditions in the K-direction are of particular importance to the solution. 
MARCHING DIRECTION I (March Direction of 1) should be chosen for flows which are 
primarily in the I-direction, and so forth. 


50 



Multiple Sweeps on Equations: In addition to choosing the solver sweep direction, you 
can control the accuracy with which each equation is solved during each iteration. This is 
accomplished by increasing the NUMBER OF SWEEPS or Number of Sweeps of an equation, 
and it is particularly useful for equations that are having difficulty converging. When you choose 
multiple sweeps for a variable, you are requesting that the equation be solved by marching across 
the domain a multiple number of times. In Fig. A. 28, this process is illustrated for the pressure 
equation, where 5 sweeps are performed by default. 

With each march across the domain, the variables at the "downstream" nodes get updated, 
thereby improving the accuracy that can be obtained during each iteration. In other words, the 
updated information obtained when a column (or row) of cells is solved gets propagated 
throughout the domain more efficiently with each additional sweep of an equation. 

Choosing the Optimum Number of Sweeps: Some guidelines are available which can 
help you choose the best number of sweeps for the various equations which are relevant to your 
problem. In problems with heat transfer, especially those with conjugate heat transfer, you should 
increase the number of sweeps on the enthalpy up to 20 or so. If you are having difficulty getting 
the pressure to converge, it is recommended that you increase the sweeps on pressure to 20 as 
well. Occasionally, you may find that one (or more) of the species equations has difficulty 
converging in some reacting flow problems. In most cases, increasing the sweeps on the 
problematic species equation to 5 or 10 is enough to make a difference. In some 3D flows, one 
velocity component may have more difficulty converging than the other two. In such cases, 
increasing the sweeps on the equation for that component to 5 or 10 is again adequate. Since one 
consequence of increasing the number of sweeps on an equation is that each iteration takes more 
time, you should request multiple sweeps for only those equations that are having difficulty 
converging. 

By default, 1 sweep is performed for each equation except the pressure-correction 
equation, for which 5 are performed. Note that these settings are not used when the equations in 
question are being solved using multigrid, as the multigrid procedure exercises its own control 
over the sweeping procedure. 

7.5.5 Block Correction 

FLUENT provides an optional block correction procedure for solution of the governing 
equations in the LGS-BLOCK-CORRECTIONS (and BLOCK-CORRECTIONS-MG) tables. 
When you are not using the multigrid option, you can add block correction on the enthalpy and 
species equations using the LGS-BLOCK-CORRECTIONS command. 

If you are using multigrid, you can add block correction on any transport equation using 
the BLOCK-CORRECTIONS-MG command. 

The block correction option accelerates the convergence of the equations by deriving 
global corrections that drive the solution toward global conservation. Block correction of the 
enthalpy equation is especially recommended for the solution of coupled convection/conduction 
problems "conjugate heat transfer" and in any problem involving regions with highly different 
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thermal conductivities. The block correction technique speeds up convergence by applying a 
quasi-one-dimensional correction to the current solution field in order to satisfy global 
conservation, with the correction vanishing as the local balance is achieved throughout the 
domain. 

ID Block Correction vs. Multigrid: One-dimensional block corrections reduce long- 
wavelength errors in the direction in which they are applied, but may introduce large short- 
wavelength errors. Therefore, application of block correction is followed by application of LGS 
to reduce the latter errors. Sufficient sweeps must be specified to this purpose, otherwise the 
effect of block correction may be adverse. If very steep gradients exist in the final solution, ID 
block correction may increase the solution effort for this reason. These disadvantages are to a 
large extent avoided in multigrid (MG) acceleration. 

One-dimensional block correction is also useful in combination with multigrid 
acceleration, when cells are highly stretched or otherwise strongly anisotropic coefficients are to 
be expected. In such cases, you can activate ID block correction with the BLOCK- 

CORRECTIONS-MG command. 

7.5.6 Multigrid Acceleration of the Line Solver 

FLUENT uses a multignd scheme to accelerate the convergence of the standard line 
solver, using additive correction to compute corrections on a series of coarse grid levels. The use 
of this multigrid scheme can greatly reduce the number of iterations and the CPU time required 
to obtain a converged solution, particularly when your FLUENT model contains a large number 
of control volumes and/or control volumes which are of large aspect ratio. This section describes 
the mathematical basis of this multigrid approach and provides details on the user inputs and 
controls over this important solver feature. 

The Line-by-Line Solver: The standard line-by-line solution technique used in FLUENT 
provides a robust and memory-efficient approach to the solution of the finite difference equations 
in FLUENT. The essence of the line-by-line solver is to update the field for a given unknown 0 
by solving simultaneous sets of equations for 0 over a subset of elements which lie on lines of 
constant I, J, or K. This approach avoids simultaneous solution of equations at all points in the 
domain, along with the large memory requirements associated with such direct matrix 
calculations. In the line-by-line approach, the matrix to be solved is a simple tri-diagonal matrix 
which can be solved via efficient Gaussian elimination. This line-by-line approach is also 
referred to as the Line-Gauss-Seidel (LGS) method. 

The Need for Multigrid Acceleration: The LGS method, in which unknown values on 
neighboring lines are left explicit during the solution process, reduces local errors with relative 
ease. That is, the effect of the solution on one line is communicated to adjacent lines relatively 
quickly. However, the line-by-line solver is less effective at reducing "long-wavelength" errors- 
errors which exist over a large number of control volumes. Thus, global corrections to the 
solution across a large number of control volumes occur slowly, over many iterations, when the 
line-by-line solver is used. This implies that performance of the line solver will deteriorate as the 
number of control volumes along any coordinate direction increases. Multigrid provides a 
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remedy for this weakness in the line solver by deriving "global" corrections which are based on a 
control volume balance over a large number of cells. Figure A. 29 shows the convergence history 
of a moderate size (36,000 nodes) 3D turbulent flow problem with and without multigrid. 

The convergence or error reduction rate of the line-by-line solver will also deteriorate 
when the finite difference coefficients in the equations to be solved are anisotropic. Anisotropic 
coefficients occur when large cell aspect ratios exist in the finite difference grid and when rapidly 
varying or anisotropic transport properties occur in the domain. One instance of anisotropic 
transport properties occurs in the energy equation when your problem includes conjugate heat 
transfer between regions of very different conductivity. Multigrid can accelerate the convergence 
of such problems by enforcing a global balance over larger regions of the grid, in effect 
"smoothing" the anisotropy that exists on the local scale. Figure A.30 illustrates the effect of 
multigrid on the solution history of a conduction heat transfer calculation in a composite 
medium, where coefficient anisotropy slows the convergence of the line solver. 

Choosing When to Use Multigrid: In most cases, multigrid (MG) should be applied to 
the pressure-correction (continuity) equation and will reduce the number of iterations required to 
converge it. Multigrid is of particular help for convergence of the continuity equation because the 
mass balance depends on reduction of long-wavelength errors. In addition, multigrid should also 
be applied to scalar equations (like enthalpy), especially when these equations are loosely 
coupled to the flow solution. In such cases, given sufficiently high underrelaxation, convergence 
of the scalar equations will primarily depend on global error reduction, and will be limited by the 
linear equation solver procedure. Again, multigrid should provide significant speedup when this 
is the case. 

Applying multigrid to momentum or "source-dominated" scalar equations, on the other 
hand, provides little if any benefit. (Source dominated scalars include the turbulence parameters 
and enthalpy or species in reacting flows.) Such equations tend to be dominated by local 
conditions, and the line-by-line solver will do well at reduction of these local errors. In fact, 
convergence may be hindered by application of multigrid to such equations as the global 
corrections may introduce significant local errors that are difficult to remove. 

7.6 Solution Options for Large Problems 

By default, FLUENT automatically updates internally stored information about your 
problem whenever you make certain changes to the problem definition. For example, when you 
make changes in the Models panel and click Apply, or when you quit from the SETUP- 1 text 
menu, FLUENT will check and update the geometry and property information that is used 
internally by the solver. (Note that this information is not saved to a file at this time- you must 
always remember to save the problem definition to a case file yourself.) 

Since FLUENT’S check and update of the geometry and properties includes verification of 
the grid, this process can be very timeconsuming for very large problems (i.e., problems with a 
very large number of cells). For such problems, you may want to enable the Prompt Before 
Geometry Calc, option in the Numerical Options panel. 
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When this option is enabled, FLUENT will prompt you before it updates the geometry 
and property information. Rather than performing this update several times during your problem 
setup, you can postpone the update until you have completed the setup. Simply answer NO in the 
text interface (or click on No in the Question dialog box) when FLUENT asks if you want to 
calculate the geometric parameters and properties. When you are ready to perform the update, 
you can either answer YES (or click on Yes) the next time you are prompted, or request the 
update yourself by selecting the Update Geometry menu item in the Solve pull-down menu. 
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Table A.l Values of 0 in common discritization schemes 


Discretization Scheme 

e 

Second-order Central Difference 

1 

Second-order Upwind 

0 

QUICK 

ffl 

Third-order Upwind 

if 


Table A.2 underrelaxation factors for SIMPLE and SIMPLEC 


Variable 

SIMPLE 

SIMPLEC 

U, U, U) 


0.7 

Swirl (u>) 

13 

0.8 

Pressure Correction 

■ffl 

1.0 

k 

0.7 

0.7 

t 

0.7 

0.7 

Enthalpy 

1.0 

1.0 

Temperature 

1.0 

1.0 

Reynolds Stresses 

0.5 

0.5 

Species 

1.0 

1.0 

Viscosity 

0.7 

0.7 

Body Forces 

1.0 

1.0 
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Fig. A.l Basic program structure. 




(b) Cartesian Coordinates 


Fig. A.2 Geometry representation in body-fitted coordinates 
vs. Cartesian coordinates. 
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(a) Approximate Geometry in Cartesian Coordinates Has a Very 
Small Impact on the Flow Prediction 



(b) Body-Fitted Coordinates Can also be Used, but Require 
More Computer Resources 


Fig. A.3 Satisfactory approximation of the geometry and flow physical coordinates. 



(a) Body-Fitted Coordinates Provide a More Accurate Descrip- 
tion of the Wall Shear Stress 


(b) Cartesian Coordinates Predict the Wall Shear Stress Inac 
curately 


Fig. A.4 Unsatisfactory approximation of the flow physics in Cartesian coordinates. 
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(a) Coarse Grid Solution 


(b) Fine Grid Solution 


Fig. A.9 Impact of the computational grid on the resolution of flow details 
in a simple driven cavity. 









Fig. A.10 Impact of the near-wall grid spacing on the shear stress 
calculation in laminar flow. 




Fig. A. 11 Undesirable skew in body-fitted grid cells. 
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Fig. A. 12 Grid lines, nodes, and control volumes in FLUENT. 
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(a) Transonic Flow in a Converging-Diverging Nozzle 


(b) Supersonic Flow Over a Bump in a 2D Channel 
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(c) Supersonic Jet Mixing in an Ambient 
Plenum 

Fig. A.16 Typical compressible flow applications modeled by FLUENMT. 



Fig. A.17 Ill-posed boundary conditions for supersonic inflow 
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Fig. A. 18 Overview of the solution process. 
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Fig. A. 19 Typical computational cell surrounding node P. 
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Fig. A.20 Non-staggered control volume storage scheme. 
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Fig. A.21 One dimensional control volume nomenclature used 
to illustrate volume integration. 
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A*U > | < AX C > | < AXp H 

Fig. A.23 Central, downwind, and upwind cell nomenclature employed 
in the higher order interpolation schemes. 
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Fig. A.24 Points and distances for viscous-weighted velocity interpolation. 
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marching direction 



Fig. A.25 Illustration of sweep direction 2 (J). 
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Fig. A.26 Illustration of sweep direction 1 (I). 




Fig. A.27 Alternating direction sweeps. 
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repeat 5 times 



times 



Fig. A.28 multiple sweeps for the pressure equation. 



Fig. A.29 Pressure-correction equation residual vs. iteration for flow 
in a 3D elbow bend. 
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